BedSet

Basic Usage

bedset.sh in=file1.bed,file2.bed out=result.bed subtract=t

Input Formats

BedSet accepts two primary input formats, automatically detected by file extension:

• BED files: Standard 3-column BED format (chrom, start, end). Coordinates are 0-based half-open intervals.
• VCF files: Automatically detected by .vcf or .vcf.gz extension. Each variant's reference span is converted to a BED interval and can be padded on each side for deletion coverage.

Coordinate System

BedSet uses 0-based half-open coordinates: an interval [start, end) has length end-start and covers no base when start==end. All input intervals are automatically sorted and merged per scaffold on load, so unsorted or self-overlapping inputs are handled correctly.

Coverage Statistics

BedSet always computes and reports base-pair coverage statistics to stderr, including:

• Total bp covered by each input
• Total bp in the union of all inputs
• Total bp in the intersection of all inputs
• For 2-input analysis: shared bp and unique bp to each input
• Output interval count and coverage

Common Use Cases

Subtract Exclusion Regions

Remove a known problematic region file from a set of genomic intervals:

bedset.sh in=regions.bed,exclusions.bed out=filtered.bed subtract=t

This command subtracts all regions in exclusions.bed from regions.bed, producing filtered.bed containing only regions not found in the exclusion set.

Create a Union of Multiple Region Sets

Merge multiple interval files into a single set of non-overlapping intervals:

bedset.sh in=set1.bed,set2.bed,set3.bed out=merged.bed union=t

This produces a single BED file with all unique regions from the three input files, with overlapping intervals merged automatically.

Find Shared Regions Between Multiple VCF Files

Identify variant positions shared across multiple VCF files with padding for deletion coverage:

bedset.sh in=variants1.vcf.gz,variants2.vcf.gz,variants3.vcf.gz out=shared.bed intersection=t pad=10

This finds regions present in all three VCF files (padded by 10 bp per side to account for deletion spans) and outputs the intersection as BED.

Build a Multiallelic-Exclusion BED from VCF

Create a BED file of positions containing multiallelic genotypes to exclude from analysis:

bedset.sh in=variants.vcf out=multiallelic.bed multiallelic=t

This extracts only variant positions where the first sample carries a multiallelic genotype (allele index ≥2), useful for building exclusion masks for downstream analyses.

Compare Coverage Between Two Region Sets

Analyze overlap statistics between two genomic interval sets:

bedset.sh in=capture.bed,exon.bed

This produces no output file but displays comprehensive statistics on bp covered by each set, the shared bp, and the unique bp to each, useful for QC of targeted sequencing designs.

Parameters

BedSet parameters are organized into several functional categories. Some parameters control the operation mode (select exactly one), others configure I/O behavior, and others handle VCF-specific processing.

Algorithm Details

Interval Processing

BedSet processes intervals using a per-scaffold coordinate sweep algorithm, as described in the source code (lines 194-258). The algorithm:

1. Loads all input intervals and stores them as parallel start/stop integer arrays per scaffold
2. Automatically sorts and merges overlapping or abutting intervals per scaffold during load
3. For each scaffold, creates boundary events (start and stop points) from all input intervals
4. Sorts events by coordinate and sweeps through them, maintaining depth counters for the first input (0 or 1) and all remaining inputs (0 to N-1)
5. At each coordinate range, applies the depth predicate for the selected operation

Set Operation Predicates

Set operations are implemented as depth predicates applied to each coordinate range (lines 266-271):

• Union: Includes coordinate ranges where depthA + depthRest ≥ 1 (covered by any input)
• Intersection: Includes ranges where depthA ≥ 1 AND depthRest ≥ N-1 (covered by all inputs)
• Subtract: Includes ranges where depthA ≥ 1 AND depthRest == 0 (covered by first input only)

VCF Variant Handling

For VCF inputs, each variant's reference span is extracted and converted to a BED interval (lines 370-401). The reference span starts at position-1 (converting from 1-based VCF to 0-based BED) and extends to position-1 + reference allele length. When pad is specified, the interval is expanded by pad bp on each side. If multiallelic=t, only variants with multiallelic genotypes in the first sample are retained.

Output Format

The output BED file is written in 3-column format (chromosome, start, end) with intervals sorted by scaffold name, then by coordinate (lines 295-310). Each interval represents a merged region from the selected set operation.

Output Statistics

BedSet always prints statistics to stderr, including the selected operation, input counts, per-input coverage information, union/intersection bp totals, and output coverage. For two-input analysis, it additionally reports shared bp and unique bp counts, with built-in self-checks validating that the sweep algebra is internally consistent (lines 160-174).

Example stderr output:

Operation:         	subtract
Inputs:            	2
Input 0 (regions.bed):
  bp covered:      	1500000
  intervals:       	450
  scaffolds:       	21
Input 1 (exclusions.bed):
  bp covered:      	250000
  intervals:       	50
  scaffolds:       	12
Union bp:          	1650000
Intersection bp:   	100000
Shared bp:         	100000
Unique to input 0: 	1400000
Unique to input 1: 	150000
Output (subtract) bp:
  bp covered:      	1400000
  intervals:       	500

Support

For questions and support:

• Email: bbushnell@lbl.gov
• Documentation: bbmap.org
• GitHub: github.com/bbushnell/BBTools