How to Convert VCF to BED: A Practical Guide

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

By the end of this guide you will be able to take a variant call file, decide deliberately which coordinate convention you need, run the conversion with a command you understand, sort and index the result, and prove that the output is correct before it reaches a downstream tool. You will also recognize the failure modes that produce an empty intersection or a silent one-base shift, because those two problems account for most of the support requests about this conversion.

Have these on hand before starting. A bgzipped and indexed VCF (bcftools, tabix). bedtools for intersection testing. The reference FASTA that matches the genome build of your VCF if you plan to normalize indels. Optionally BEDOPS for its vcf2bed utility and Python 3 with pysam if you prefer a scripted approach. All commands below use short flags so they paste cleanly into a terminal.

VCF is described in the literature as the lingua franca of next-generation sequencing, which is a fair summary of its position in genomics workflows [1]. BED is the lingua franca of genomic intervals. The two formats describe overlapping information with different coordinate systems, and the conversion is simple only if you make the coordinate decision explicitly.

What VCF and BED Actually Store

VCF: one record per variant

A VCF file is tab-delimited text. Metadata lines begin with two hash characters. The column header line begins with a single hash followed by CHROM. After the header, each line is one variant record with eight fixed fields, an optional FORMAT field, and one field per sample.

Three properties of VCF cause most conversion mistakes. POS is 1-based. POS refers to the first base of the reference allele, not to the center of the variant. And a single record can carry several alternate alleles at once, written G,T in the ALT field, so one line does not always mean one variant.

VCF also carries sample-level genotypes, phasing, and depth. Those are the parts that do not survive a conversion to BED, and they are the reason you should keep the original VCF as the source of truth. VCF is a rich format that has been used as the input side of many downstream conversions, including clinical conversions into HL7 FHIR Genomics reports where each VCF row becomes a structured variant record [1]. The BED format has no equivalent capacity.

BED: one record per interval

BED is also tab-delimited. It requires three columns and permits up to twelve.

The single most important rule in this entire article is that BED counts from zero and excludes its end coordinate. A variant at 1-based position 100 occupies the interval from 0-based 99 up to but not including 100. It is one base wide. Writing chr1 100 100 produces a zero-width interval that most tools silently ignore, and writing chr1 100 101 shifts the variant one base to the right.

Where FASTQ and FASTA fit

These four formats sit at different points in a sequencing pipeline, which is why confusing them causes trouble.

FASTQ stores raw reads in four lines per record: a header beginning with @, the nucleotide sequence, a separator line beginning with +, and a per-base quality string in Phred+33 encoding. FASTA stores sequences in two lines per record: a header beginning with > and the sequence itself, with no quality information at all. FASTA is the usual container for reference genomes and protein sequences. FASTQ is the usual container for instrument output.

VCF and BED come later. A VCF records positions where a sample differs from the reference. A BED records arbitrary intervals with no sequence and no alleles. Neither can be reconstructed from the other. Reading a FASTQ file as if it were FASTA silently folds quality strings into sequence records, which is exactly the class of error that well-built format libraries are designed to prevent by validating structure rather than trusting the extension [2].

The Coordinate Decision That Determines Everything

One base, worked two ways

Take this VCF record:

chr1	100	rs1234	A	G	60	PASS	DP=30

Read it as follows. On chromosome 1, at 1-based position 100, the reference base is A and one alternate base is G. The variant covers exactly one reference base.

The equivalent BED interval is:

chr1	99	100

The arithmetic is start = POS - 1 and end = start + length(REF), which gives 99 + 1 = 100.

Now take a deletion:

chr1	100	.	ATG	A	50	PASS	.

The reference allele is three bases long. The variant covers 1-based positions 100, 101, and 102. The equivalent BED interval is:

chr1	99	102

The arithmetic is start = 100 - 1 = 99 and end = 99 + 3 = 102. This is the reference span, and it is the convention that bcftools reports with its %END placeholder.

Four interval conventions and when each applies

There is no universally correct answer to "how wide is a variant in BED." There are four conventions in common use, and each is right for a particular question.

The point interval is the simplest and the most common in tutorials. It is fine when your VCF contains only single nucleotide variants. It becomes misleading for indels, because a three-base deletion represented as a one-base interval will miss annotation features that overlap the deleted bases.

The reference span is the safest general-purpose choice. It matches what bedtools assumes when it reads a VCF directly, and it matches the %END placeholder in bcftools query.

The alternate span matters when you care about inserted sequence rather than deleted sequence. For an insertion such as REF=A, ALT=ATTT, the reference span is one base while the alternate span is four bases.

Whichever convention you choose, write it down in the methods section of your analysis. A BED file carries no record of how it was built.

Multiallelic records

A record like this:

chr1	100	.	A	G,T	45	PASS	.

describes two alternate alleles at one position. Converting it directly yields one interval for the reference allele, which does not distinguish the two alternates. If you need one interval per allele, split first:

bcftools norm -m -any -f reference.fa -Oz -o split.vcf.gz input.vcf.gz
bcftools index -t split.vcf.gz

The -m -any option splits multiallelic sites into biallelic records. After splitting, each record has a single ALT and the interval arithmetic becomes unambiguous.

Chromosome naming and record order

Two BED files can look identical and share zero intervals because one uses chr1 and the other uses 1. This produces no error message. bedtools intersect returns an empty file and exits successfully, which is the most expensive kind of failure.

Check the naming in the VCF header before converting:

bcftools view -h input.vcf.gz | grep '^##contig' | head -5

The output lists the contig names the caller used. Match those names in every BED file you plan to intersect against. If you must rename, do it once, early, and keep the map file.

Record order matters too. Most interval tools assume coordinate-sorted input. Sorting with plain sort on the whole line is wrong because it sorts the start column as text. Use the numeric key form:

sort -k1,1 -k2,2n variants.bed > variants.sorted.bed

How to Convert VCF to BED Step by Step

Step 1: Inventory the input

Never convert a file you have not looked at.

bcftools view -h input.vcf.gz | head -30
bcftools view -H input.vcf.gz | head -3
bcftools index -n input.vcf.gz

The first command prints the metadata header. The second prints the first three records so you can see whether IDs are present and whether ALT contains commas. The third prints the total number of records, which you will use later as a checksum.

A typical header fragment looks like this:

##fileformat=VCFv4.2
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant">

Seeing the ##contig lines confirms both the naming scheme and the presence of length metadata. Seeing an ##INFO=<ID=END line tells you that symbolic alleles in this file will carry an explicit end position.

Step 2: Normalize and split

Normalization left-aligns indels against the reference and splits multiallelic sites. Both steps reduce ambiguity in the conversion.

bcftools norm -f reference.fa -m -any -Oz -o normalized.vcf.gz input.vcf.gz
bcftools index -t normalized.vcf.gz

The -f option supplies the reference FASTA. Without it, norm can still split records but cannot left-align indels or verify that the REF allele matches the genome. The -Oz option writes compressed VCF output and -o names it. The -t flag on bcftools index builds a tabix index rather than a CSI index, which is what most downstream tools expect for VCF.

The output is silent on success. Check the record count changed if splitting occurred:

bcftools index -n normalized.vcf.gz

Step 3: Run the conversion

The fastest and most reliable conversion for a compressed VCF uses bcftools query with a format string:

bcftools query -f '%CHROM\t%POS0\t%END\n' normalized.vcf.gz > variants.bed

Three placeholders do the work. %CHROM prints the contig name. %POS0 prints the position converted to 0-based, which is exactly the BED start. %END prints the end of the reference allele, which is the 0-based exclusive end. The tab escapes inside the format string keep the output column-delimited.

The first lines of the output look like this:

chr1	99	100
chr1	149	151
chr1	299	300
chr2	499	499

The first record is a single-base variant at 1-based 100. The second is a two-base reference span starting at 1-based 150, so it covers 150 and 151. The last record has identical start and end, which indicates a zero-width or malformed entry and deserves investigation before you proceed.

For a point-interval conversion instead, replace %END with %POS:

bcftools query -f '%CHROM\t%POS0\t%POS\n' normalized.vcf.gz > variants.point.bed

Step 4: Add name and score columns

A three-column BED is often enough, but annotation tools frequently want a label. Add the ID:

bcftools query -f '%CHROM\t%POS0\t%END\t%ID\n' normalized.vcf.gz > variants.bed

Where the ID is a dot, substitute a constructed name from the coordinates:

bcftools query -f '%CHROM\t%POS0\t%END\t%ID\t%QUAL\n' normalized.vcf.gz \
  | awk 'BEGIN{OFS="\t"} { if ($4 == ".") $4 = $1 ":" ($2 + 1); print }' \
  > variants.named.bed

The $2 + 1 rebuilds the 1-based position so the name reads like chr1:100, which is what most genome browsers display. The QUAL field lands in the score column, which is legal only if your qualities fall between 0 and 1000. Check that before you rely on it.

Step 5: Sort and index

sort -k1,1 -k2,2n variants.named.bed > variants.sorted.bed
bgzip -c variants.sorted.bed > variants.sorted.bed.gz
tabix -p bed variants.sorted.bed.gz

The -p bed preset tells tabix that the file uses BED coordinates and that the first three columns are chromosome, start, and end. The result is variants.sorted.bed.gz and variants.sorted.bed.gz.tbi. The index enables region queries later.

Step 6: Validate before use

Four checks catch nearly every error.

wc -l variants.sorted.bed
bcftools index -n normalized.vcf.gz

These two numbers should match if you emitted one interval per record. A mismatch after splitting means some records were dropped, and a mismatch before splitting means you collapsed multiallelic sites.

awk 'BEGIN{OFS="\t"} $3 <= $2 {print "zero or negative width:", $0}' variants.sorted.bed | head

Any output from this command means an interval has no width. Zero-width records are almost always a coordinate bug.

bcftools view -H -r chr1:100-100 normalized.vcf.gz
grep -P '^chr1\t99\t100\t' variants.sorted.bed

The first command prints the VCF record at 1-based position 100. The second prints the matching BED lines. Confirm they refer to the same variant. Note the shift in the region syntax: bcftools view -r takes 1-based inclusive coordinates while BED stores 0-based half-open coordinates.

bedtools intersect -a variants.sorted.bed -b genes.bed -wa | wc -l

A non-zero overlap count against a known gene model confirms that chromosome naming and sorting are compatible. Zero overlaps against a gene set you know should overlap is the signal to check the chr prefix.

Step 7: Use the result

bedtools intersect -a variants.sorted.bed -b genes.bed -wa -wb > variants.in.genes.bed

The -wa flag writes the original variant interval. The -wb flag appends the overlapping gene interval. The output is a joined table you can pass to downstream counting or annotation.

Command Cookbook

The awk approach is a useful fallback when bcftools is unavailable, but it assumes the file is uncompressed and that POS is the second column, which is true for VCF and false for VCF-like tables exported from other software. It also ignores the reference allele length, so it only produces point intervals.

The BEDOPS route is fast and memory-efficient for large files. Piping through sort-bed guarantees the coordinate sort that BEDOPS tools expect.

The pysam route gives you programmatic control:

import pysam

vcf = pysam.VariantFile("normalized.vcf.gz")
with open("variants.bed", "w") as out:
    for rec in vcf:
        label = rec.id if rec.id else "."
        out.write(f"{rec.chrom}\t{rec.start}\t{rec.stop}\t{label}\n")

In pysam, rec.start is the 0-based start of the reference allele and rec.stop is the exclusive end. Those are precisely the BED start and end, so no manual arithmetic is needed. This is the advantage of using a library that models the format rather than treating it as text. Libraries built around validated format models exist across bioinformatics, from quantum chemistry to proteomics, and they exist precisely because ad hoc text munging breaks at scale [2][3].

Working With Large VCFs

bcftools query streams records and does not load the file into memory, so a multi-gigabyte VCF converts with modest RAM. For repeated work, convert once to BCF, the binary counterpart of VCF, and query that:

bcftools view -Ob -o input.bcf input.vcf.gz
bcftools index input.bcf
bcftools query -f '%CHROM\t%POS0\t%END\n' input.bcf > variants.bed

Random access requires bgzip compression plus a tabix or CSI index. A plain gzip file cannot be queried by region. Writing a BED without sorting first also wastes time later, because virtually every interval tool either requires sorted input or re-sorts internally.

When the conversion is one step in a longer pipeline, wrapping it in a workflow engine pays off. Workflow libraries provide job automation, reproducibility, and audit traces for every output file, which matters when a coordinate convention decision needs to be documented alongside the results [4]. Dedicated toolkits also tend to outperform generic converters as data volume grows, as benchmarks for Hi-C format conversion have shown [5]. The same principle applies here: for a one-off conversion, bcftools query is enough, and for a production pipeline, an explicit, tested step is worth the setup.

Common Mistakes and Limitations

Forgetting the minus one. The single most common error. VCF POS is 1-based and BED start is 0-based. Every conversion command in this guide handles it, but hand-written awk that prints $2 in the start column does not.

Writing POS in both start and end columns. chr1 100 100 has zero width. Most tools drop it without complaint. The correct point interval is chr1 99 100.

Using POS as the end after subtracting one from the start. chr1 99 101 is two bases wide and shifts every variant one position to the right. Compute the end as start + length(REF), not as POS + 1.

Mismatched chromosome naming. chr1 and 1 do not match. bedtools intersect returns an empty file with exit status zero. There is no warning. Check the ##contig header first.

Unsorted input. bedtools with its sorted mode reports that the file has an out of order record and names the offending line. Fix with sort -k1,1 -k2,2n, not with a plain sort.

Spaces instead of tabs. A BED file delimited by spaces fails with a malformed BED entry error and an illegal field count. Set OFS="\t" in awk or use bcftools query, which emits tabs.

Collapsed multiallelic sites. One record with two ALTs becomes one interval. If your downstream analysis counts alleles rather than sites, split with bcftools norm -m -any first.

Symbolic alleles producing one-base intervals. A record with ALT=<DEL> has a one-character REF allele, so the reference-span calculation yields a single base. Use the INFO/END field instead:

bcftools query -f '%CHROM\t%POS0\t%INFO/END\n' normalized.vcf.gz > sv.bed

Records lacking an END value will fail or emit a placeholder. Filter them out first or handle them separately.

Indexing the wrong file. Running tabix on a plain gzip file or on an uncompressed BED produces an htslib open failure. Compress with bgzip, not gzip.

Confusing BED coordinates with region syntax. BED is 0-based half-open. bcftools view -r and samtools view take 1-based inclusive regions. A BED interval chr1 99 100 is chr1:100-100 in region syntax, not chr1:99-100.

Expecting a round trip. BED cannot store REF, ALT, genotype, phase, quality, or filter status. Converting BED back to VCF produces a syntactically valid file with no allele information, which is worse than useless for most analyses because it looks like data. Keep the VCF.

Assuming the conversion is lossless for structural variants. Breakend notation, symbolic alleles, and multi-allelic structural records do not have a clean interval representation. Treat these separately and document what you did.

Formatting errors that survive to the next step. Conversion tools that skip validation push the failure downstream, where it is harder to diagnose. GUI wrappers around converters reduce typing errors but do not eliminate the need to inspect output, a point made explicit by tools that convert aligned FASTA into population-genetics input formats [6]. Standardized, human-readable formats with explicit schemas are easier to validate than ad hoc text, which is the argument made for newer formats in mass spectrometry and beyond [7]. Applying the same discipline here means checking counts, widths, and a handful of spot coordinates every time.

One more limitation worth stating plainly. Interoperability between formats is an active engineering problem, not a solved one, and even bidirectional converters between well-documented formats require careful handling of proprietary extensions [8]. VCF to BED is simpler than most such problems because both formats are open and well specified, but the coordinate mismatch is a real semantic gap. Scripted conversions that move data between formats need explicit tests, and that is true whether the target is a sequence submission checklist [9] or an interval file [10].

Frequently Asked Questions

What is the difference between VCF POS and BED start?

VCF POS is 1-based and points at the first base of the reference allele. BED chromStart is 0-based. For a variant at POS 100, the BED start is 99.

Does converting VCF to BED lose genotype information?

Yes, completely. BED has no fields for genotype, allele depth, phase, or filter status. Only coordinates and optional labels survive, so keep the VCF as the authoritative record.

Should I subtract 1 from POS when making a BED file?

Yes, always. Subtract 1 for the start and compute the end as start plus the length of the reference allele. Commands like bcftools query -f '%POS0' do the subtraction for you.

How do I convert only SNPs to BED?

Filter first, then convert. Use bcftools view -v snps input.vcf.gz to keep only single nucleotide variants, then run the standard query. Filtering by variant type also removes the indel ambiguity from the coordinate decision.

Why does bedtools intersect return zero overlaps with my VCF-derived BED?

The usual cause is a chromosome naming mismatch, such as chr1 in one file and 1 in the other. The second most common cause is inconsistent sort order. Neither produces an error message.

Can I convert BED back to VCF?

Only syntactically. BED carries no REF or ALT alleles, so a reconstructed VCF has no variant identity. Round-tripping through BED is a data-loss operation, not a format conversion.

What is the difference between FASTQ and FASTA?

FASTQ stores reads in four lines with a per-base quality string. FASTA stores sequences in two lines with no quality data. FASTQ is instrument output and FASTA is typically a reference or an assembled sequence.

How do I convert a VCF with structural variants to BED?

Use the INFO/END field rather than the reference allele length. Query with %INFO/END and handle records that lack it separately, since symbolic alleles do not have a reference span that reflects their true size.

Sources

1. vcf2fhir: a utility to convert VCF files into HL7 FHIR format for genomics-EHR integration.
2. IOData: A python library for reading, writing, and converting computational chemistry file formats and generating input files.
3. ThermoRawFileParser: Modular, Scalable, and Cross-Platform RAW File Conversion.
4. SciPipe: A workflow library for agile development of complex and dynamic bioinformatics pipelines.
5. hictk: blazing fast toolkit to work with .hic and .cool files.
6. Fasta2Structure: a user-friendly tool for converting multiple aligned FASTA files to STRUCTURE format.
7. HXMS: a standardized file format for HX-MS data.
8. cd2sbgnml: bidirectional conversion between CellDesigner and SBGN formats.
9. EMBL2checklists: A Python package to facilitate the user-friendly submission of plant and fungal DNA barcoding sequences to ENA.
10. anndataR improves interoperability between R and Python in single-cell transcriptomics.