%prog refine breakpoints.bed gaps.bed

    Find gaps within or near breakpoint region.

    For breakpoint regions with no gaps, there are two options:
    - Break in the middle of the region
    - Break at the closest gap (--closest)
160664
def refine(args):
    
    p = OptionParser(refine.__doc__)
    p.add_option("--closest", default=False, action="store_true",
                 help="In case of no gaps, use closest [default: %default]")
    opts, args = p.parse_args(args)

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

    breakpointsbed, gapsbed = args
    ncols = len(open(breakpointsbed).next().split())
    logging.debug("File {0} contains {1} columns.".format(breakpointsbed, ncols))
    cmd = "intersectBed -wao -a {0} -b {1}".format(breakpointsbed, gapsbed)

    pf = "{0}.{1}".format(breakpointsbed.split(".")[0], gapsbed.split(".")[0])
    ingapsbed = pf + ".bed"
    sh(cmd, outfile=ingapsbed)

    fp = open(ingapsbed)
    data = [x.split() for x in fp]

    nogapsbed = pf + ".nogaps.bed"
    largestgapsbed = pf + ".largestgaps.bed"
    nogapsfw = open(nogapsbed, "w")
    largestgapsfw = open(largestgapsbed, "w")
    for b, gaps in groupby(data, key=lambda x: x[:ncols]):
        gaps = list(gaps)
        gap = gaps[0]
        if len(gaps) == 1 and gap[-1] == "0":
            assert gap[-3] == "."
            print("\t".join(b), file=nogapsfw)
            continue

        gaps = [(int(x[-1]), x) for x in gaps]
        maxgap = max(gaps)[1]
        print("\t".join(maxgap), file=largestgapsfw)

    nogapsfw.close()
    largestgapsfw.close()
    beds = [largestgapsbed]
    toclean = [nogapsbed, largestgapsbed]

    if opts.closest:
        closestgapsbed = pf + ".closestgaps.bed"
        cmd = "closestBed -a {0} -b {1} -d".format(nogapsbed, gapsbed)
        sh(cmd, outfile=closestgapsbed)
        beds += [closestgapsbed]
        toclean += [closestgapsbed]
    else:
        pointbed = pf + ".point.bed"
        pbed = Bed()
        bed = Bed(nogapsbed)
        for b in bed:
            pos = (b.start + b.end) / 2
            b.start, b.end = pos, pos
            pbed.append(b)
        pbed.print_to_file(pointbed)
        beds += [pointbed]
        toclean += [pointbed]

    refinedbed = pf + ".refined.bed"
    FileMerger(beds, outfile=refinedbed).merge()

    # Clean-up
    FileShredder(toclean)

    return refinedbed
1243475
def location(args):
    
    from jcvi.formats.bed import BedLine
    from jcvi.graphics.histogram import stem_leaf_plot

    p = OptionParser(location.__doc__)
    p.add_option("--dist", default=100, type="int",
                 help="Distance cutoff to call 5` and 3` [default: %default]")
    opts, args = p.parse_args(args)

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

    bedfile, fastafile = args
    dist = opts.dist
    sizes = Sizes(fastafile).mapping
    fp = open(bedfile)
    fiveprime = threeprime = total = 0
    percentages = []
    for row in fp:
        b = BedLine(row)
        pos = b.start
        size = sizes[b.seqid]
        if pos < dist:
            fiveprime += 1
        if size - pos < dist:
            threeprime += 1
        total += 1
        percentages.append(100 * pos / size)

    m = "Five prime (within {0}bp of start codon): {1}\n".format(dist, fiveprime)
    m += "Three prime (within {0}bp of stop codon): {1}\n".format(dist, threeprime)
    m += "Total: {0}".format(total)
    print(m, file=sys.stderr)

    bins = 10
    title = "Locations within the gene [0=Five-prime, 100=Three-prime]"
    stem_leaf_plot(percentages, 0, 100, bins, title=title)
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
714267
def setBreakpoints( self, breakpoints ):
        
        self.clearBreakpoints()
        for breakpoint in breakpoints:
            self.addBreakpoint(breakpoint)
348343
def paste(args):
    
    from jcvi.formats.bed import uniq

    p = OptionParser(paste.__doc__)
    p.add_option("--maxsize", default=300000, type="int",
            help="Maximum size of patchers to be replaced [default: %default]")
    p.add_option("--prefix", help="Prefix of the new object [default: %default]")
    p.set_rclip(rclip=1)
    opts, args = p.parse_args(args)

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

    pbed, blastfile, bbfasta = args
    maxsize = opts.maxsize  # Max DNA size to replace gap
    order = Bed(pbed).order

    beforebed, afterbed = blast_to_twobeds(blastfile, order, log=True,
                                           rclip=opts.rclip, maxsize=maxsize,
                                           flipbeds=True)
    beforebed = uniq([beforebed])

    afbed = Bed(beforebed)
    bfbed = Bed(afterbed)

    shuffle_twobeds(afbed, bfbed, bbfasta, prefix=opts.prefix)
1317171
def overlap(args):
    
    from jcvi.apps.console import green

    p = OptionParser(overlap.__doc__)
    p.add_option("--maxerr", default=2, type="int", help="Maximum error rate")
    p.add_option("--canvas", default=100, type="int", help="Canvas size")
    opts, args = p.parse_args(args)

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

    bestcontains, iid = args
    canvas = opts.canvas

    bestcontainscache = bestcontains + ".cache"
    if need_update(bestcontains, bestcontainscache):
        fp = open(bestcontains)
        fw = open(bestcontainscache, "w")
        exclude = set()
        for row in fp:
            if row[0] == '#':
                continue
            j = int(row.split()[0])
            exclude.add(j)
        dump(exclude, fw)
        fw.close()

    exclude = load(open(bestcontainscache))
    logging.debug("A total of {0} reads to exclude".format(len(exclude)))

    cmd = "overlapStore -d ../asm.ovlStore -b {0} -e {0}".format(iid)
    cmd += " -E {0}".format(opts.maxerr)
    frags = []
    for row in popen(cmd):
        r = OverlapLine(row)
        if r.bid in exclude:
            continue
        frags.append(r)

    # Also include to query fragment
    frags.append(OverlapLine("{0} {0} N 0 0 0 0".format(iid)))
    frags.sort(key=lambda x: x.ahang)

    # Determine size of the query fragment
    cmd = "gatekeeper -b {0} -e {0}".format(iid)
    cmd += " -tabular -dumpfragments ../asm.gkpStore"
    fp = popen(cmd)
    row = next(fp)
    size = int(fp.next().split()[-1])

    # Determine size of canvas
    xmin = min(x.ahang for x in frags)
    xmax = max(x.bhang for x in frags)
    xsize = -xmin + size + xmax
    ratio = xsize / canvas

    fw = sys.stdout
    for f in frags:
        fsize = -f.ahang + size + f.bhang
        a = (f.ahang - xmin) / ratio
        b = fsize / ratio
        t = '-' * b
        if f.orientation == 'N':
            t = t[:-1] + '>'
        else:
            t = '<' + t[1:]
        if f.ahang == 0 and f.bhang == 0:
            t = green(t)
        c = canvas - a - b
        fw.write(' ' * a)
        fw.write(t)
        fw.write(' ' * c)
        print("{0} ({1})".format(str(f.bid).rjust(10), f.erate_adj), file=fw)
649387
def frombed(args):
    
    p = OptionParser(frombed.__doc__)
    p.add_option("--gapsize", default=100, type="int",
                 help="Insert gaps of size [default: %default]")
    opts, args = p.parse_args(args)

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

    bedfile, = args
    gapsize = opts.gapsize
    agpfile = bedfile.replace(".bed", ".agp")
    fw = open(agpfile, "w")

    bed = Bed(bedfile, sorted=False)
    for object, beds in groupby(bed, key=lambda x: x.accn):
        beds = list(beds)
        for i, b in enumerate(beds):
            if gapsize and i != 0:
                print("\t".join(str(x) for x in \
                    (object, 0, 0, 0, "U", \
                     gapsize, "scaffold", "yes", "map")), file=fw)

            print("\t".join(str(x) for x in \
                    (object, 0, 0, 0, "W", \
                     b.seqid, b.start, b.end, b.strand)), file=fw)

    fw.close()

    # Reindex
    return reindex([agpfile, "--inplace"])
1161177
def bed(args):
    
    p = OptionParser(bed.__doc__)
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    pslfile, = args
    fw = must_open(opts.outfile, "w")

    psl = Psl(pslfile)
    for p in psl:
        print(p.bed12line, file=fw)
79401
def insert(args):
    
    from jcvi.formats.agp import mask, bed
    from jcvi.formats.sizes import agp

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

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

    candidates, gapsbed, chrfasta, unplacedfasta = args
    refinedbed = refine([candidates, gapsbed])
    sizes = Sizes(unplacedfasta).mapping
    cbed = Bed(candidates)
    corder = cbed.order
    gbed = Bed(gapsbed)
    gorder = gbed.order

    gpbed = Bed()
    gappositions = {}  # (chr, start, end) => gapid

    fp = open(refinedbed)
    gap_to_scf = defaultdict(list)
    seen = set()
    for row in fp:
        atoms = row.split()
        if len(atoms) <= 6:
            continue
        unplaced = atoms[3]
        strand = atoms[5]
        gapid = atoms[9]
        if gapid not in seen:
            seen.add(gapid)
            gi, gb = gorder[gapid]
            gpbed.append(gb)
            gappositions[(gb.seqid, gb.start, gb.end)] = gapid
        gap_to_scf[gapid].append((unplaced, strand))

    gpbedfile = "candidate.gaps.bed"
    gpbed.print_to_file(gpbedfile, sorted=True)

    agpfile = agp([chrfasta])
    maskedagpfile = mask([agpfile, gpbedfile])
    maskedbedfile = maskedagpfile.rsplit(".", 1)[0] + ".bed"
    bed([maskedagpfile, "--outfile={0}".format(maskedbedfile)])

    mbed = Bed(maskedbedfile)
    finalbed = Bed()
    for b in mbed:
        sid = b.seqid
        key = (sid, b.start, b.end)
        if key not in gappositions:
            finalbed.add("{0}\n".format(b))
            continue

        gapid = gappositions[key]
        scfs = gap_to_scf[gapid]

        # For scaffolds placed in the same gap, sort according to positions
        scfs.sort(key=lambda x: corder[x[0]][1].start + corder[x[0]][1].end)
        for scf, strand in scfs:
            size = sizes[scf]
            finalbed.add("\t".join(str(x) for x in \
                    (scf, 0, size, sid, 1000, strand)))

    finalbedfile = "final.bed"
    finalbed.print_to_file(finalbedfile)

    # Clean-up
    toclean = [gpbedfile, agpfile, maskedagpfile, maskedbedfile]
    FileShredder(toclean)
217588
def bed(args):
    
    p = OptionParser(bed.__doc__)
    p.add_option("--switch", default=False, action="store_true",
                 help="Switch reference and aligned map elements [default: %default]")
    opts, args = p.parse_args(args)

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

    mapout, = args
    pf = mapout.split(".")[0]
    mapbed = pf + ".bed"
    bm = BinMap(mapout)
    bm.print_to_bed(mapbed, switch=opts.switch)

    return mapbed
420528
public static void main(final String[] args) {
        Switch about = new Switch("a", "about", "display about message");
        Switch help = new Switch("h", "help", "display help message");
        FileArgument bedFile = new FileArgument("o", "bed-file", "output BED file, default stdout", false);
        IntegerArgument n = new IntegerArgument("n", "n", "number of BED records to generate, default " + DEFAULT_N, false);
        IntegerArgument size = new IntegerArgument("s", "size", "chromosome size, default " + DEFAULT_SIZE, false);
        StringArgument chrom = new StringArgument("c", "chrom", "chromosome name, default " + DEFAULT_CHROM, false);
        DoubleArgument meanLength = new DoubleArgument("l", "mean-length", "mean length, default " + DEFAULT_MEAN_LENGTH, false);
        DoubleArgument lengthVariation = new DoubleArgument("v", "length-variation", "length variation, default " + DEFAULT_LENGTH_VARIATION, false);
        IntegerArgument seed = new IntegerArgument("z", "seed", "random number seed, default relates to current time", false);

        ArgumentList arguments = new ArgumentList(about, help, bedFile, n, size, chrom, meanLength, lengthVariation, seed);
        CommandLine commandLine = new CommandLine(args);

        GenerateBed generateBed = null;
        try
        {
            CommandLineParser.parse(commandLine, arguments);
            if (about.wasFound()) {
                About.about(System.out);
                System.exit(0);
            }
            if (help.wasFound()) {
                Usage.usage(USAGE, null, commandLine, arguments, System.out);
                System.exit(0);
            }

            RandomGenerator random = seed.wasFound() ? new MersenneTwister(seed.getValue()) : new MersenneTwister();
            double lv = Math.max(NO_VARIATION, lengthVariation.getValue(DEFAULT_LENGTH_VARIATION));
            RealDistribution length = new NormalDistribution(random, meanLength.getValue(DEFAULT_MEAN_LENGTH), lv, NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);

            generateBed = new GenerateBed(bedFile.getValue(), n.getValue(DEFAULT_N), size.getValue(DEFAULT_SIZE), chrom.getValue(DEFAULT_CHROM), random, length);
        }
        catch (CommandLineParseException e) {
            Usage.usage(USAGE, e, commandLine, arguments, System.err);
            System.exit(-1);
        }
        try {
            System.exit(generateBed.call());
        }
        catch (Exception e) {
            e.printStackTrace();
            System.exit(1);
        }
    }
416497
def tidy(args):
    
    p = OptionParser(tidy.__doc__)
    p.add_option("--gapsize", dest="gapsize", default=0, type="int",
            help="Set all gaps to the same size [default: %default]")
    p.add_option("--minlen", dest="minlen", default=100, type="int",
            help="Minimum component size [default: %default]")

    opts, args = p.parse_args(args)

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

    fastafile, = args
    gapsize = opts.gapsize
    minlen = opts.minlen

    tidyfastafile = fastafile.rsplit(".", 1)[0] + ".tidy.fasta"
    fw = must_open(tidyfastafile, "w")

    removed = normalized = 0
    fasta = Fasta(fastafile, lazy=True)
    for name, rec in fasta.iteritems_ordered():
        rec.seq = rec.seq.upper()
        if minlen:
            removed += remove_small_components(rec, minlen)
        trim_terminal_Ns(rec)
        if gapsize:
            normalized += normalize_gaps(rec, gapsize)

        if len(rec) == 0:
            logging.debug("Drop seq {0}".format(rec.id))
            continue
        SeqIO.write([rec], fw, "fasta")

    # Print statistics
    if removed:
        logging.debug("Total discarded bases: {0}".format(removed))
    if normalized:
        logging.debug("Gaps normalized: {0}".format(normalized))

    logging.debug("Tidy FASTA written to `{0}`.".format(tidyfastafile))
    fw.close()

    return tidyfastafile
569835
private static String[] cutOffAtHelp(String[] args) {
		List<String> argsToKeep = new ArrayList<>();
		for (String arg : args) {
			argsToKeep.add(arg);
			if (arg.endsWith(CliCommandLineProducer.OPTION_HELP_POSTFIX)) {
				break;
			}
		}
		
		if (argsToKeep.isEmpty()) {
			argsToKeep.add("--help");
		}
		
		return argsToKeep.toArray(new String[]{});
	}
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)
381115
def at_breakpoint(breakpoint)
      @pry ||= Pry.new

      output.puts bold("\n  Breakpoint #{breakpoint.id}. ") + n_hits(breakpoint)

      expr = breakpoint.expr
      return unless expr

      output.puts bold("Condition: ") + expr
    end
754773
def fixboundaries(args):
    
    p = OptionParser(fixboundaries.__doc__)
    p.add_option("--type", default="gene", type="str",
                 help="Feature type for which to adjust boundaries")
    p.add_option("--child_ftype", default="mRNA", type="str",
                 help="Child featuretype(s) to use for identifying boundaries")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    gffile, = args
    gffdb = make_index(gffile)

    fw = must_open(opts.outfile, "w")
    for f in gffdb.all_features(order_by=('seqid', 'start')):
        if f.featuretype == opts.type:
            child_coords = []
            for cftype in opts.child_ftype.split(","):
                for c in gffdb.children(f, featuretype=cftype, order_by=('start')):
                    child_coords.append((c.start, c.stop))
            f.start, f.stop = range_minmax(child_coords)

        print(f, file=fw)

    fw.close()
405804
def tiling(args):
    
    p = OptionParser(tiling.__doc__)
    p.add_option("--overlap", default=3000, type="int",
                 help="Minimum amount of overlaps required")
    p.set_verbose()
    opts, args = p.parse_args(args)

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

    bedfile, = args
    ov = opts.overlap

    bed = Bed(bedfile)
    inf = len(bed)
    selected = Bed()
    for seqid, sbed in bed.sub_beds():
        g = Grouper()
        current = sbed[0]
        # Partition connected features
        for a in sbed:
            g.join(a)
            # requires a real overlap
            if a.start < current.end - ov:
                g.join(a, current)
            if a.end > current.end:
                current = a

        # Process per partition
        for gbed in g:
            end = max(x.end for x in gbed)
            gbed.sort(key=lambda x: (x.start, -x.end))
            entries = len(gbed)
            counts = [inf] * entries
            counts[0] = 1
            traceback = [-1] * entries
            for i, a in enumerate(gbed):
                for j in xrange(i + 1, entries):
                    b = gbed[j]
                    if b.start >= a.end - ov:
                        break
                    # Two ranges overlap!
                    if counts[i] + 1 < counts[j]:
                        counts[j] = counts[i] + 1
                        traceback[j] = i
            endi = [i for i, a in enumerate(gbed) if a.end == end]
            last = min((traceback[i], i) for i in endi)[1]
            chain = []
            while last != -1:
                chain.append(last)
                last = traceback[last]
            chain = chain[::-1]
            selected.extend([gbed[x] for x in chain])

            if opts.verbose:
                print(counts)
                print(traceback)
                print(chain)
                print("\n".join(str(x) for x in gbed))
                print("*" * 30)
                print("\n".join(str(gbed[x]) for x in chain))
                print()

    tilingbedfile = bedfile.rsplit(".", 1)[0] + ".tiling.bed"
    selected.print_to_file(filename=tilingbedfile, sorted=True)
    logging.debug("A total of {} tiling features written to `{}`"\
                    .format(len(selected), tilingbedfile))