How to Sort and Merge BED Files: Best Practices

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-10

By the end of this guide you will be able to take any pile of BED files, put them into a reproducible sort order, merge overlapping intervals without losing the columns you care about, and prove with a couple of cheap checks that the output is correct. You need a Unix-like shell, GNU coreutils (the sort that ships with most Linux systems and with Homebrew coreutils on macOS), and BEDTools. SAMtools is optional unless you are starting from BAM files. BEDOPS is optional but worth installing if your intervals run into the hundreds of millions.

The commands below are complete and runnable. Every one of them has been written to avoid hidden dependencies, and every expected output is shown so you can compare against your own run.

What BED Files Are and Why Order Matters

The BED coordinate system

BED is a plain-text format that describes genomic intervals on a reference assembly. Its coordinates are zero-based and half-open, meaning the first base of a chromosome is position 0 and an interval covering bases 1 through 100 is written as 0 100. The end coordinate is exclusive, so the length of an interval is simply end - start. A chr1 0 100 record covers 100 bases.

This convention differs from GFF, GTF, and VCF, all of which are one-based. A GFF feature written 1 100 covers the same 100 bases as the BED record 0 100. If you convert formats without subtracting one from the start, every interval shifts by a single base. That shift is small enough to survive visual inspection and large enough to break overlap analysis at boundaries. The Ensembl GFF3 specification is explicit that GFF coordinates are one-based and inclusive, so the subtraction is not optional.

BED3 through BED12

BED is defined in a series of column counts, and the count determines how the file is interpreted.

• BED3 contains chromosome, start, and end. It is the minimum required by bedtools merge.
• BED4 adds a name in column 4.
• BED5 adds a score in column 5.
• BED6 adds strand in column 6, using +, -, or ..
• BED12 adds thickStart, thickEnd, itemRgb, blockCount, blockSizes, and blockStarts, which together describe a spliced feature such as a transcript.

The block fields in BED12 are relative: blockStarts are offsets from the interval start, not absolute chromosome coordinates. bedtools merge reads the first three columns and ignores everything else unless you explicitly ask for more with -c and -o. A BED12 file that you merge will come back as BED3, and any exon structure is gone. That behavior surprises people who expect merge to preserve transcripts.

Header, track, and browser lines

UCSC-style BED files often begin with lines that start with #, track, or browser. These are display directives, not data. BEDTools will either reject them or misparse them depending on the subcommand, and GNU sort will happily sort them into the middle of your data. Strip them before sorting.

BED is tab-delimited. A file that uses runs of spaces instead of tabs will parse incorrectly in strict tools and produce column-count errors downstream. Check this at ingest.

Step 1: Inspect and Clean the Input

Run these checks before you sort anything. They take seconds and catch the majority of problems that would otherwise surface as confusing downstream errors.

1. Confirm the field count is consistent.

awk -F'\t' '{print NF}' input.bed | sort -n | uniq -c

A clean BED6 file produces a single line such as:

  18422 6

More than one count means some lines have a different number of columns. That is usually a header line, a blank line, or a line with spaces instead of tabs.

1. Count comment and track lines.

grep -c -E '^(#|track|browser)' input.bed

1. Check for Windows line endings.

grep -c $'\r' input.bed

Any nonzero count means the file has carriage returns. Convert with tr -d '\r' before sorting. A trailing \r on the last column ruins numeric sort keys and makes the end coordinate unparseable.

1. Check for spaces used as delimiters.

grep -c ' ' input.bed

1. Check coordinate sanity.

awk -F'\t' '$2 < 0 || $3 <= $2' input.bed | head

This prints intervals with a negative start or an end that is not greater than the start. Both are invalid BED.

1. Inspect chromosome names.

cut -f1 input.bed | sort -u | head -30

Look for a mixture of chr1 and 1, or chrM and MT. Those two naming systems are both common, and mixing them produces zero overlaps rather than an error.

1. Produce a clean working copy.

grep -v -E '^(#|track|browser)' input.bed | tr -d '\r' > clean.bed

Keep the original untouched. Every transformation after this point should be reproducible from clean.bed.

Step 2: How to Sort BED Files

The canonical sort command

This is the command to memorize:

LC_ALL=C sort -k1,1 -k2,2n -k3,3n clean.bed > sorted.bed

Each piece does specific work.

• LC_ALL=C forces byte-order collation. Without it, sort behavior depends on the locale, and the same file can sort differently on two machines. Reproducibility across runs and across colleagues depends on this variable.
• -k1,1 sorts on the first field only, the chromosome.
• -k2,2n sorts on the second field, the start, using numeric comparison.
• -k3,3n breaks ties on the end coordinate so the output is deterministic when two intervals share a start.
• -s requests a stable sort, which preserves the original relative order of records with identical keys. Add it if downstream logic depends on input order.

bedtools merge requires input sorted by chromosome and then by start position within each chromosome. It does not require any particular chromosome order. It does require that all records for one chromosome are contiguous, which the command above guarantees.

Why the n matters

Dropping the n from -k2,2n is the single most common sorting error, and it produces a file that looks plausible while being wrong.

Consider this unsorted file:

chr2	50	150	peak_3	20	-
chr1	100	250	peak_1	50	+
chr1	1000	1200	peak_4	10	+
chr1	200	400	peak_2	80	+

Sorting with a text comparison on the start column gives:

chr1	100	250	peak_1	50	+
chr1	1000	1200	peak_4	10	+
chr1	200	400	peak_2	80	+
chr2	50	150	peak_3	20	-

The string 1000 sorts before the string 200 because 1 comes before 2 character by character. Feeding that file to bedtools merge produces:

Error: Sorted input specified, but the file sorted_bad.bed has the following out of order record
chr1	200	400	peak_2	80	+

The fix is the numeric suffix. The correct output is:

chr1	100	250	peak_1	50	+
chr1	200	400	peak_2	80	+
chr1	1000	1200	peak_4	10	+
chr2	50	150	peak_3	20	-

If you ever see the "out of order record" error, the cause is almost always a missing n, a file that was never sorted, or two files concatenated without re-sorting.

Chromosome order: lexicographic versus natural

With LC_ALL=C, chromosomes sort as bytes: chr1, chr10, chr11, chr2. That order is stable and perfectly acceptable to BEDTools, as long as every file in a pipeline uses the same one. Problems appear when one input was sorted lexicographically and another was sorted naturally. The concatenated file then has chromosome blocks interleaved, and merge fails or, worse, produces fragmented output.

Three options give you explicit control.

LC_ALL=C sort -k1,1 -k2,2n -k3,3n clean.bed > sorted_lex.bed
LC_ALL=C sort -k1,1V -k2,2n -k3,3n clean.bed > sorted_natural.bed
bedtools sort -faidx reference.fa.fai -i clean.bed > sorted_ref.bed

The first is the default recommendation. The second uses version sort, which is available in GNU coreutils and gives chr1, chr2, chr3, chr10. macOS ships BSD sort, which does not support -V, so install GNU coreutils through Homebrew and call gsort if you need natural order there.

The third option is the most defensible when a reference assembly is available. bedtools sort -faidx reads the order of contigs from the .fai index produced by samtools faidx and sorts the BED accordingly. Every file sorted against the same index ends up in exactly the same chromosome order, which removes a whole class of pipeline bugs.

BEDTools sort, GNU sort, or sort-bed

All three produce valid input for bedtools merge. They differ in speed, control, and validation.

bedtools sort is the simplest. It sorts by chromosome and start, accepts a genome file with -g or a FASTA index with -faidx, and is convenient inside existing BEDTools pipelines.

GNU sort is the fastest option once you tune memory and temporary storage. It gives you locale control, tie-breaking, and stable sorting. For anything above a few million intervals, this is the tool to use.

sort-bed from BEDOPS is a specialized BED sorter that validates the format as it sorts and defines a chromosome ordering designed for genomic data. It is a good choice for very large inputs and for pipelines that later use bedops, merge.

Sorting very large files

Sorting is the expensive step in this workflow. Merging is a streaming operation and costs almost nothing by comparison. Give the sorter memory and a fast scratch disk.

LC_ALL=C sort -k1,1 -k2,2n -k3,3n -S 25% -T /scratch/tmp big.bed > big.sorted.bed

-S 25% tells GNU sort to use up to 25 percent of physical memory before spilling to disk. -T /scratch/tmp puts the spill files on local disk rather than a network mount. On a machine with 64 GB of RAM, this uses roughly 16 GB and usually avoids disk spilling for interval files in the tens of millions of records.

Write the sorted file to disk rather than piping it directly into merge if you intend to reuse it. Re-sorting the same file three times in one project is a waste of hours.

Sorting with strand

For stranded data, add strand as a sort key so that records with the same coordinates but opposite strands are grouped predictably.

LC_ALL=C sort -k1,1 -k2,2n -k6,6 -k3,3n stranded.bed > stranded.sorted.bed

This assumes column 6 exists and contains strand. Sorting this way costs nothing and makes downstream -s merges easier to reason about.

Step 3: How to Merge BED Files

What merge actually does

bedtools merge walks a sorted BED file from start to finish and collapses intervals that meet the merge criteria into single records. The output is BED3 by default, containing chromosome, start, and end.

bedtools merge -i sorted.bed > merged.bed

The semantics are simple. Records are considered for merging when they are on the same chromosome and within the distance threshold. The merged interval extends from the smallest start to the largest end in the group.

The distance parameter

-d controls how close two intervals must be to merge. It is the parameter people misunderstand most often.

• -d 0 is the default. It merges intervals that overlap and also merges intervals that are book-ended, meaning the end of one exactly equals the start of the next.
• -d -1 requires at least one base of true overlap. Book-ended intervals stay separate.
• -d 1000 merges intervals separated by up to 1000 bases, which is useful for collapsing nearby peaks into regions but dangerous if you want to preserve distinct loci.

Given this sorted input:

chr1	100	250	peak_1	50	+
chr1	200	400	peak_2	80	+
chr1	1000	1200	peak_4	10	+
chr2	50	150	peak_3	20	-

The default merge produces:

chr1	100	400
chr1	1000	1200
chr2	50	150

The first two intervals overlap and collapse. The gap between 400 and 1000 is 600 bases, well beyond the default threshold, so the third interval stays separate. Running with -d 1000 instead produces:

chr1	100	1200
chr2	50	150

The distinction matters. A merge that is too aggressive hides real structure, and a merge that is too conservative leaves you with the same redundancies you started with. Decide the biological question first, then choose -d.

Keeping and summarizing columns

Merging discards everything beyond the first three columns unless you ask for it. The -c flag names the columns to retain or summarize, and -o names the operation to apply to each.

bedtools merge -i sorted.bed -c 4,5 -o count,max > summarized.bed

Expected output for the example above:

chr1	100	400	2	80
chr1	1000	1200	1	10
chr2	50	150	1	20

The number of entries in -c must match the number of entries in -o. A mismatch produces a column and operation count error, and the message names both lists so you can correct it.

Common operations include count, distinct, collapse, sum, mean, median, min, max, first, and last. distinct keeps unique values separated by commas, while collapse keeps every value including duplicates. Change the separator with -delim, for example -delim "|".

bedtools merge -i sorted.bed -c 4 -o distinct -delim "|" > names.bed

This is the standard way to attach sample names to merged regions. It is also the standard way people accidentally create enormous column values when a merge collapses thousands of intervals. Check the longest field in the output before you load it into a spreadsheet.

Stranded merging

By default, merge ignores strand. A + interval and a - interval that overlap become a single unstranded record. For RNA-seq clusters, antisense transcripts, or any strand-specific assay, that is incorrect.

bedtools merge -i stranded.sorted.bed -s -c 4 -o count > stranded.merged.bed

The -s flag requires column 6 and merges only features on the same strand. The output gains a strand column. The -S flag does the opposite and merges only features on opposite strands, which is rarely what you want but exists for specific analyses.

Merging several files at once

The principle is always the same: combine, sort the combined stream, then merge. Sorting each file individually and concatenating the results is not enough, because the combined file is not sorted.

cat sample_a.bed sample_b.bed sample_c.bed \
  | LC_ALL=C sort -k1,1 -k2,2n -k3,3n \
  | bedtools merge -i - > union.bed

The -i - argument tells BEDTools to read from standard input.

If you need to know which samples contributed to each merged region, use multiinter instead.

bedtools multiinter -i sample_a.bed sample_b.bed sample_c.bed -header > presence.bed

The output contains chromosome, start, end, the number of files overlapping the interval, a comma-separated list of those files, and one indicator column per file. The -header flag adds column names. This is the right tool for a Venn-diagram-style question, and merge is the right tool when you want a single non-redundant set of regions.

For a fully merged union with per-sample annotation, run multiinter and then merge the first three columns again. Do not try to do both jobs with one command.

Worked Example, Start to Finish

This walkthrough uses two small files so you can verify every step.

Start with two unsorted inputs.

sample_a.bed:

chr1	100	250	peak_a1	50	+
chr1	1000	1200	peak_a4	10	+
chr2	50	150	peak_a3	20	-

sample_b.bed:

chr1	200	400	peak_b1	80	+
chr1	1050	1300	peak_b2	30	+
chr2	500	700	peak_b3	15	-

Step 1. Combine and sort.

cat sample_a.bed sample_b.bed \
  | LC_ALL=C sort -k1,1 -k2,2n -k3,3n \
  > combined.sorted.bed

cat combined.sorted.bed

chr1	100	250	peak_a1	50	+
chr1	200	400	peak_b1	80	+
chr1	1000	1200	peak_a4	10	+
chr1	1050	1300	peak_b2	30	+
chr2	50	150	peak_a3	20	-
chr2	500	700	peak_b3	15	-

Step 2. Merge and keep a count.

bedtools merge -i combined.sorted.bed -c 4 -o count > combined.merged.bed

cat combined.merged.bed

chr1	100	400	2
chr1	1000	1300	2
chr2	50	150	1
chr2	500	700	1

Step 3. Merge with a 2 kb window to collapse nearby loci on each chromosome.

bedtools merge -i combined.sorted.bed -d 2000 -c 4 -o count > combined.wide.bed

cat combined.wide.bed

chr1	100	1300	4
chr2	50	700	2

Step 4. Confirm the merged file is stable under a second merge.

bedtools merge -i combined.merged.bed | wc -l
bedtools merge -i combined.merged.bed -c 4 -o count | wc -l

Both counts must equal the original line count of combined.merged.bed, which is 4. If they do not, the first merge was not applied to sorted input.

Key Options at a Glance

Verifying the Result

A merge that runs without errors can still be wrong. Four checks catch nearly every failure.

1. Re-merge check. Merging an already merged file must not change the line count or the coordinates.

bedtools merge -i merged.bed | diff - merged.bed

An empty diff means the file is idempotent under merge. Any output means the original merge was incomplete.

1. Coordinate check. No interval should have a negative start or an end less than or equal to its start.

awk -F'\t' '$2 < 0 || $3 <= $2' merged.bed | wc -l

The expected count is zero.

1. Coverage check. Sum the interval lengths before and after.

awk -F'\t' '{s += $3 - $2} END {print s}' clean.bed
awk -F'\t' '{s += $3 - $2} END {print s}' merged.bed

The merged total is the true number of covered bases. The input total is an overcount whenever the input contains overlaps. The difference tells you how much redundancy the merge removed.

1. Overlap check against the source. Every input interval should be contained in some merged interval.

bedtools intersect -a clean.bed -b merged.bed -v | head

The -v flag reports input intervals with no overlap in the merged file. An empty result, or one that contains only intervals you expected to lose, confirms the merge covered the input.

For stranded data, repeat the checks with -s and a matching sort. For very large files, use bedtools jaccard -a merged.bed -b reference.bed to quantify similarity to a known region set rather than checking by eye.

Common Mistakes and Limitations

Sorting with -k2,2 instead of -k2,2n. This is the classic error. It produces a file that looks sorted to a human and fails inside merge. Always use the numeric suffix on coordinate keys.

Sorting by start without grouping by chromosome. A sort on column 2 alone produces a file where chromosomes are interleaved. Merge then treats the whole file as one unsorted stream.

Mixing chr1 and 1. Human-readable references use one or the other. UCSC uses the chr prefix and Ensembl does not, and mitochondrial contigs appear as chrM in one and MT in the other. Overlaps between mismatched names are zero, not an error. Normalize at ingest with sed 's/^/chr/' or sed 's/^chr//', and handle the mitochondrial rename explicitly.

Leaving header lines in place. track and browser lines break parsers. Strip them first.

Spaces instead of tabs. BED is tab-delimited. Files exported from spreadsheets often use spaces or mixed whitespace. Rebuild the file with a strict tab delimiter.

Carriage returns. A file that was edited on Windows carries \r at the end of every line. The end coordinate then fails numeric parsing, and the sort key becomes a string.

Merging stranded data without -s. Overlapping plus and minus features collapse into one record. For strand-specific work this is a silent biological error, not a formatting one.

Merging across feature classes. Exons, genes, and peaks are different objects. Merging them together produces intervals that do not correspond to any real annotation. Split by feature type first.

Expecting merge to preserve BED12 structure. Merge returns BED3 by default. If you need exon-level merging, expand the blocks first with bedtools bed12tobed6 -i transcripts.bed, merge the resulting BED6 blocks, and rebuild transcripts afterward if needed.

Choosing -d carelessly. A 10 kb window can collapse genuinely distinct regulatory loci. Pick the distance from the biology, not from the desire to reduce the line count.

Concatenating sorted files without re-sorting. Sorting each file separately does not sort the concatenation. Always pipe the combined stream through a single sort.

Assuming merge counts duplicates. Merge collapses geometry, not records. If two identical intervals exist, -c 4 -o count reports 2. Use sort -u on BED3 before merging if you want identical intervals removed.

Forgetting the genome file for complement. bedtools complement needs a two-column file of chromosome names and lengths. Generate it from a FASTA index.

samtools faidx reference.fa
cut -f1,2 reference.fa.fai > genome.txt

Extra columns. BED permits additional columns beyond 12, and merge will tolerate them, but other subcommands may reject them. Trim to the columns you need before running a mixed pipeline.

The limitations are equally worth stating. bedtools merge has no concept of biological meaning. It merges coordinates, and it will merge a promoter with a distal enhancer if the distance threshold allows it. It cannot merge across chromosome boundaries. It cannot infer strand from a file without a strand column. It does not fill gaps, and it does not know about assembly gaps unless you run complement against a genome file. Every one of those behaviors is a deliberate design choice, and each one has to be handled by the analyst.

Where BED Fits Among FASTQ, FASTA, BAM, and GFF

FASTQ versus FASTA

These two formats are frequently confused because they both hold sequence, and neither one is an interval format.

FASTA stores sequence records with a > header line followed by the sequence, usually wrapped across multiple lines. It carries no quality information. Reference genomes, protein sequences, and transcript sets are distributed as FASTA.

FASTQ stores sequence and per-base quality in four lines per record: an @ header, the sequence, a + separator, and a quality string of the same length as the sequence. The quality characters are Phred scores encoded with an offset of 33 in the current Illumina and Sanger convention. FASTQ is the raw output of a sequencer, and the quality string is what allows base calling errors to be modeled.

The practical difference is what you can do with each. A FASTQ file tells you how confident the sequencer was in each base. A FASTA file tells you what the sequence is. A BED file tells you where on a reference something happened.

Converting between them is routine.

awk 'NR % 4 == 1 {print ">" substr($0, 2)} NR % 4 == 2 {print}' reads.fastq > reads.fasta

The seqtk utility does the same job with more options.

seqtk seq -A reads.fastq > reads.fasta

From BAM to BED

Alignment files are the usual source of BED intervals. Coordinate-sort the BAM first, then convert.

samtools sort -@ 8 -o aligned.sorted.bam aligned.bam
bedtools bamtobed -i aligned.sorted.bam > reads.bed

samtools sort orders records by reference and position, which is the same ordering bedtools merge expects. The resulting BED therefore arrives already sorted. If you filter the BAM afterward with samtools view, the order is preserved. If you combine BAMs with samtools merge, the inputs must be sorted and the output is sorted, so the converted BED still needs no extra sort. Verify with a re-merge check anyway.

Extracting sequence with a BED and a FASTA

bedtools getfasta pulls the sequence under each interval from a reference FASTA.

samtools faidx reference.fa
bedtools getfasta -fi reference.fa -bed merged.bed -fo merged.fa

The -s flag respects strand and reverse-complements features on the minus strand. -name uses column 4 as the output record name, and -tab writes a tab-delimited table instead of FASTA. The FASTA index is mandatory. Without it, getfasta fails with a message about the missing .fai file.

Converting GFF or GTF to BED

GFF is one-based, so the start coordinate needs the subtraction described earlier.

awk 'BEGIN{OFS="\t"} $0 !~ /^#/ {print $1, $4-1, $5, $3, $6, $7}' annotations.gff > annotations.bed

That produces BED6 from a standard nine-column GFF3. Check the strand column and the feature types before merging, because a GFF typically mixes genes, transcripts, and exons in one file.

Practical Tips for Production Pipelines

Normalize chromosome names once, at ingest. Pick chr-prefixed or unprefixed, convert every input, and keep the mapping in the pipeline. Do not normalize in three different scripts.

Keep the reference index and the genome file with the data. A .fai file and a two-column genome file let you sort in reference order and run complement without guessing contig lengths.

Re-sort after any transformation that changes coordinates. Intersections, subtractions, slflanking operations, and shifts can all disturb order. A sort before every merge is cheap insurance.

Stream instead of writing intermediates. For large jobs, pipe sort directly into merge.

LC_ALL=C sort -k1,1 -k2,2n -k3,3n input.bed | bedtools merge -i - > merged.bed

Log the line count at every stage. A sudden drop of 90 percent between two steps means a filter or a merge did something you did not intend.

Use strict shell settings in scripts. set -euo pipefail stops a pipeline at the first failure instead of letting a truncated file flow downstream.

Compress and index the final product. A bgzipped, tabix-indexed BED supports random access without decompressing the whole file.

bgzip merged.sorted.bed
tabix -p bed merged.sorted.bed.gz
tabix merged.sorted.bed.gz chr1:1000000-1100000

tabix requires the file to be sorted and bgzipped. The -p bed preset sets the column and coordinate conventions for you.

Version the reference assembly. A BED file without an assembly label is ambiguous. Put hg38 or GRCh38 in the filename and in the pipeline metadata.

Test on a small slice first. Run the full pipeline on one chromosome before launching it on the whole genome. Errors in sort order and chromosome naming show up immediately at small scale.

Frequently Asked Questions

Do BED files have to be sorted before merging?

Yes. bedtools merge requires input sorted by chromosome and then by start position. Unsorted input produces an out-of-order record error or, worse, fragmented output that looks valid.

What is the difference between sort -k2,2 and sort -k2,2n?

The n suffix makes the sort numeric. Without it, sort compares coordinates as text, so 1000 is placed before 200. The numeric form is required for any coordinate column.

Does bedtools merge join book-ended intervals?

Yes, by default. With -d 0, intervals that overlap and intervals where one end exactly equals the next start are both merged. Use -d -1 to require at least one base of true overlap.

How do I merge only intervals that truly overlap?

Run bedtools merge -d -1 -i sorted.bed. The negative distance requires a real overlap and leaves book-ended intervals as separate records.

How do I keep the name or score column after merging?

Use -c to list the columns and -o to list the operations, for example bedtools merge -i sorted.bed -c 4,5 -o distinct,max. The two lists must have the same length.

How do I merge BED files from different samples?

Concatenate them, sort the combined stream, and merge. cat a.bed b.bed | LC_ALL=C sort -k1,1 -k2,2n -k3,3n | bedtools merge -i - produces a single non-redundant interval set. Use bedtools multiinter if you need per-sample presence information.

What is the difference between FASTQ and FASTA?

FASTQ carries per-base quality scores in a four-line record, while FASTA carries sequence only in a two-line record. FASTQ is raw sequencing output, and FASTA is used for reference genomes and other sequences without quality data.

Can I merge BED files that use different chromosome naming?

Not meaningfully. A file using chr1 and a file using 1 will show zero overlaps. Rename one convention to match the other before sorting and merging, and handle chrM versus MT separately.