How to Convert FASTA to BED: A Step-by-Step Guide

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

By the end of this guide you will be able to turn a FASTA file into a valid, validated BED file, and you will know exactly which of the four legitimate conversion scenarios you are in. You will be able to build a chromosome-size BED from an assembly, recover intervals from coordinate-bearing FASTA headers, encode a gapped alignment as BED12 blocks, and check every output line against the source FASTA before you hand the file to bedtools, IGV, or a pipeline.

You need four things on hand: the FASTA file itself, the name of the reference assembly it came from, a decision about where the coordinates are supposed to come from, and a small toolchain. Install samtools, bedtools, seqkit, and Python with pyfaidx or Biopython. Any recent version works. Distribution-packaged builds through Conda and BioContainers are the standard way to keep converter versions pinned across machines and pipeline runs [1].

What FASTA and BED Actually Store

FASTA holds letters, not positions

A FASTA file is plain text. Each record starts with a header line beginning with > and continues with sequence lines that are usually wrapped at 60, 70, or 80 characters. The record identifier is the first whitespace-delimited token of the header. Everything after that first token is free-text description.

>chr1
ACGTACGTACGTACGTACGT
>chr2
TTTTGGGGCCCCAAAA

That file contains two sequences of 20 bp and 16 bp. It contains no information about where those sequences sit in a genome, how they relate to each other, or what strand they are on. A finished human genome FASTA is roughly 3.1 billion characters of this kind of text, and every one of those characters is unanchored to a coordinate system.

BED holds positions, not letters

BED is the UCSC interval format. Each line describes one feature on one chromosome, with three required tab-separated columns: chrom, chromStart, chromEnd. Coordinates are 0-based and half-open. The start is included in the interval, the end is excluded.

chr1	0	1

That single line is the first base of chr1. The line chr1 0 20 covers exactly 20 bases, from position 1 through position 20 in 1-based counting. Optional columns add a name, a score from 0 to 1000, a strand, thick-start and thick-end for drawing, an RGB color, and block fields. BED12 stores blockCount, blockSizes, and blockStarts, which is how a spliced transcript or a gapped alignment is represented.

BED is strictly line-oriented. There is no line wrapping in BED. One feature occupies exactly one line. This is the single most important structural difference between the two formats: FASTA wraps freely, BED cannot wrap at all.

The gap that defines the conversion problem

Because FASTA has no coordinates and BED has nothing but coordinates, "convert FASTA to BED" is never a pure format translation. It is always one of four operations:

If none of those four applies, you do not have a FASTA-to-BED problem. You may have a BED-to-FASTA problem, which is what bedtools getfasta solves, or a FASTQ problem, which is different again.

Conversion between related formats is routinely harder than a split on whitespace, and dedicated toolkits exist because general-purpose parsing is not enough. Interoperability between two formats that both describe the same underlying data still forces purpose-built converters when each format carries its own extensions and conventions [2], and format-specific quirks are the usual reason a naive parser fails [3]. The standard architecture for this work is a library that parses a file into an in-memory record, then writes the target format from that record, keeping reading and writing as separate steps [4]. Plain-text formats survive because they are readable and easy to inspect, a design trade-off that newer genomics-inspired formats still copy [5].

What You Need Installed

A minimal environment for everything in this guide:

mamba create -n fasta2bed -c conda-forge -c bioconda samtools bedtools seqkit pyfaidx biopython
mamba activate fasta2bed

Run each tool with no arguments to read its version banner, then record the versions next to your output files. Version drift is the most common reason a conversion that worked last year produces a different line count this year.

The tools and what each is for:

• samtools faidx builds a .fai index from a FASTA and gives you exact sequence lengths.
• awk turns the .fai or the headers into BED columns.
• seqkit fx2tab exports names and lengths without an index.
• bedtools getfasta runs the reverse direction and is your round-trip check.
• pyfaidx and Biopython read FASTA records in Python. Biopython in particular is the standard way to pull an identifier and a sequence out of a FASTA record separately [6].

Step-by-Step: Genome FASTA to a Whole-Chromosome BED

This is the most common request. You have genome.fa and you want a BED file listing every chromosome from position 0 to its full length.

Step 1: Inspect the FASTA before touching it

Three commands answer three questions: how many records, how the file is wrapped, and whether the identifiers are unique.

grep -c '^>' genome.fa
head -n 4 genome.fa
awk '/^>/ {print $1}' genome.fa | sort | uniq -d

For the toy file above, the first command prints 2. The second prints the first header and the first sequence line, which tells you the wrapping width. The third prints nothing when identifiers are unique and prints the repeated identifier when they are not.

Pay attention to the header style. >chr1 is clean. >chr1 Homo sapiens chromosome 1 is also usable, because samtools faidx uses only the first token. >chr1|human and >1 both work but will not match annotation files that use a different naming convention.

Step 2: Build the FASTA index

samtools faidx genome.fa

This creates genome.fa.fai, a tab-delimited file with five columns per record.

For the toy file:

chr1	20	6	20	21
chr2	16	33	16	17

The last three columns depend on header length and wrapping, so your numbers will differ. Column 2 is the one you care about, and it is exact.

samtools faidx will refuse to index a FASTA with irregular line lengths inside a single record, reporting an error such as [E::fai_build_core] Different line length in sequence 'chr1'. Re-wrap the file to a fixed width and rerun.

Step 3: Convert the index to BED

awk 'BEGIN{OFS="\t"} {print $1, 0, $2}' genome.fa.fai > genome.bed
cat genome.bed

Expected output:

chr1	0	20
chr2	0	16

The arithmetic is worth stating explicitly. A chromosome of length L is represented as [0, L). The number of bases covered is L - 0, which is L. On a real GRCh38 assembly the same command yields a first line of chr1 0 248956422, and the sum of the third column across all records equals the total assembly length.

Run the whole thing as one line if you prefer:

samtools faidx genome.fa && awk 'BEGIN{OFS="\t"} {print $1, 0, $2}' genome.fa.fai > genome.bed

Step 4: Validate the BED

Four checks, each a one-liner. A correct BED passes all of them.

awk 'NF != 3 {print "line " NR ": expected 3 columns, found " NF}' genome.bed

awk '$2 >= $3 {print "line " NR ": start >= end: " $0}' genome.bed

awk '
NR == FNR { LEN[$1] = $2; next }
{ if (!($1 in LEN)) print "unknown chrom:", $0 }
{ if ($1 in LEN && $3 > LEN[$1]) print "end past chrom end:", $0 }
' genome.fa.fai genome.bed

awk '{s += $3 - $2} END {print s " bases in BED"}' genome.bed
awk '{s += $2} END {print s " bases in FASTA index"}' genome.fa.fai

Empty output from the first three means no structural problems. The fourth check prints two totals that must be identical. If the BED total is smaller, you have dropped records. If it is larger, at least one end coordinate is wrong.

Step 5: Confirm with a round trip

bedtools getfasta turns BED back into FASTA. If your BED is correct, the regenerated sequence lengths must match the originals.

bedtools getfasta -fi genome.fa -bed genome.bed -fo roundtrip.fa
grep -v '^>' roundtrip.fa | tr -d '\n' | wc -c

The count should equal the sum from Step 4. This is the strongest single check because it exercises both the FASTA and the BED through an independent tool.

Step 6: Use the result

bedtools makewindows -g genome.fa.fai -w 1000000 > windows.bed

makewindows accepts the two-column genome file that the first two .fai columns already provide. For windowing work this is more direct than building a whole-chromosome BED first. For coverage normalization, scaffold definitions, or anything that wants one interval per chromosome, keep genome.bed.

Alternative paths without samtools

If you do not want an index, seqkit reads the FASTA directly:

seqkit fx2tab -n -l genome.fa | awk 'BEGIN{OFS="\t"} {print $1, 0, $2}' > genome.bed

In Python, pyfaidx builds the index on first access:

from pyfaidx import Fasta

genome = Fasta("genome.fa")
with open("genome.bed", "w") as out:
    for record in genome:
        out.write(f"{record.name}\t0\t{len(record)}\n")

Biopython does the same without any index, at the cost of reading the whole file:

from Bio import SeqIO

with open("genome.fa") as handle, open("genome.bed", "w") as out:
    for record in SeqIO.parse(handle, "fasta"):
        out.write(f"{record.id}\t0\t{len(record.seq)}\n")

Note the difference between record.id and record.description. The first is the identifier, the second is the full header line. Writing the description into the BED name column produces spaces in a tab-delimited file, which breaks downstream parsers.

Step-by-Step: Coordinate-Encoded Headers to BED

Regions extracted from a genome often come back with the coordinates baked into the header. bedtools getfasta produces headers of the form >chr1:1001-2000. If your FASTA looks like this, you can recover the intervals.

>chr1:1001-2000
ACGT...
>chr2:5001-5500
TTTT...

Step 1: Decide whether the numbers are 1-based or 0-based

This is the decision that causes most wrong BED files. Two conventions exist:

• 1-based, inclusive: chr1:1001-2000 is 1000 bases, the first of which is base 1001.
• 0-based, half-open: this is BED's own convention and appears in headers written by tools that pass BED coordinates straight through.

You cannot tell from the text alone. Verify it. Fetch the same region with samtools faidx, which uses 1-based inclusive region strings, and compare the first few bases to the FASTA record.

samtools faidx genome.fa chr1:1001-1005
head -n 2 regions.fa

If the five bases match the beginning of the record, the header is 1-based and you subtract 1 from the start. If they match the record shifted by one, the header is already BED-style and you subtract nothing.

Step 2: Parse the header into three columns

awk 'BEGIN{OFS="\t"}
     /^>/ {
       header = substr($0, 2)
       n = split(header, parts, ":")
       split(parts[n], se, "-")
       print parts[4], se[4] - 1, se[3]
     }' regions.fa > regions.bed

Expected output for the example above:

chr1	1000	2000
chr2	5000	5500

The conversion formula is simple and worth memorizing:

The end coordinate does not change in either convention, because an inclusive 1-based end and an exclusive 0-based end describe the same final base. Only the start needs adjusting.

Step 3: Handle names that contain a colon or a dash

seq-1:100-200 breaks naive splitting if the name itself contains a dash. A regex anchored at the end of the header is safer:

import re
import sys

pattern = re.compile(r"^(?P<name>.+?):(?P<start>\d+)-(?P<end>\d+)$")

for line in open("regions.fa"):
    if line.startswith(">"):
        header = line[1:].strip()
        match = pattern.match(header)
        if match is None:
            sys.exit(f"Cannot parse coordinates from header: {header}")
        start = int(match["start"])
        end = int(match["end"])
        print(f"{match['name']}\t{start - 1}\t{end}")

The lazy .+? lets the name contain colons and dashes while the anchors force the final colon and final dash to be the coordinate separators.

Step 4: Add strand if the header carries it

Headers such as >chr1:1001-2000(-) carry orientation. Append the strand as column 6:

awk 'BEGIN{OFS="\t"}
     /^>/ {
       header = substr($0, 2)
       strand = "+"
       if (header ~ /\(-$\)/) strand = "-"
       gsub(/\([-+]\)$/, "", header)
       n = split(header, parts, ":")
       split(parts[n], se, "-")
       print parts[4], se[4] - 1, se[3], ".", 0, strand
     }' regions.fa > regions_bed6.bed

That produces valid BED6 with a placeholder name and score. Replace the . and 0 if you have real values. A BED score must be an integer between 0 and 1000. Fractional scores will be rejected or silently truncated by different tools.

A Single Python Script for Both Cases

Most real work mixes both scenarios in one batch. This script handles whole-record intervals and header-encoded intervals with a mode flag.

#!/usr/bin/env python3
"""Convert a FASTA file to BED in one of two modes."""

import re
import sys
from pyfaidx import Fasta

HEADER_COORD = re.compile(r"^(?P<name>.+?):(?P<start>\d+)-(?P<end>\d+)$")

def whole_records(path):
    genome = Fasta(path)
    for record in genome:
        yield record.name, 0, len(record)

def header_regions(path):
    genome = Fasta(path)
    for record in genome:
        match = HEADER_COORD.match(record.name)
        if match is None:
            sys.exit(f"Header {record.name!r} has no start-end coordinates")
        start = int(match["start"])
        end = int(match["end"])
        if start < 1:
            sys.exit(f"Header {record.name!r} has a 0-based start coordinate")
        if end < start:
            sys.exit(f"Header {record.name!r} has end before start")
        yield match["name"], start - 1, end

def main():
    if len(sys.argv) != 3:
        sys.exit("usage: fasta_to_bed.py {sizes|regions} input.fa")
    mode, path = sys.argv[4], sys.argv[3]
    stream = whole_records(path) if mode == "sizes" else header_regions(path)
    for name, start, end in stream:
        print(f"{name}\t{start}\t{end}")

if __name__ == "__main__":
    main()

Run it and redirect:

python fasta_to_bed.py sizes genome.fa > genome.bed
python fasta_to_bed.py regions regions.fa > regions.bed

The start < 1 guard is deliberate. It catches the common case where someone pastes 0-based coordinates into a script that subtracts one, producing an invalid negative start. Failing loudly at that point is much cheaper than debugging a negative coordinate three steps downstream.

How to Verify the BED Is Right

Run these checks in order. Each one catches a different class of error.

1. Line count matches the FASTA record count: grep -c '^>' genome.fa against wc -l < genome.bed.
2. Every line has the expected number of columns: the NF check above.
3. Every start is less than its end: the $2 >= $3 check above.
4. Every end is at or below the chromosome length: the .fai join check above.
5. The total covered bases equal the expected total.
6. The file is tab-delimited with no header row: head -n 1 genome.bed | cat -A shows ^I between fields if they are tabs.
7. The file ends with a newline: tail -c 1 genome.bed | od -c shows \n.
8. A round trip through bedtools getfasta regenerates the expected sequence lengths.
9. Spot checks: samtools faidx genome.fa chr1:1-10 compared against the first line of the regenerated FASTA.

If your BED will feed a tool that requires sorted input, sort it and test for sortedness before the tool complains:

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

A pipeline that converts FASTA to BED should keep the conversion step scripted and logged, so that every BED file traces back to the exact FASTA and parameter set that produced it. Workflow libraries that record a full audit trace for each output file exist specifically to make that traceability automatic [7].

Common Mistakes and Limitations

Treating the BED start as 1-based

BED starts at 0. If you write chr1 1 20 intending the first 20 bases, you have described bases 2 through 20, which is 19 bases. The line chr1 0 20 is the correct description. Every off-by-one error in BED traces back to this.

Using an inclusive end coordinate

BED ends are exclusive. The interval chr1 0 20 covers 20 bases, not 21. If you copy an end coordinate from a 1-based inclusive source without checking, every interval gains one base and any length comparison fails by exactly one per feature.

Writing a header row

A BED file has no header. The line chrom start end at the top of the file will cause bedtools to stop with an error reporting an invalid record and echoing the offending line. UCSC tolerates track and browser lines, and many tools tolerate # comments, but a plain column-name header is not part of the format.

Space-delimited columns

BED requires tabs. Spaces sometimes parse, sometimes do not, and fail silently when a name field contains a space. cat -A is the fastest way to see what separates your fields.

Assuming FASTA headers are unique

Duplicate identifiers cause samtools faidx to build an index with ambiguous lookups and cause region extraction to return the wrong record. awk '/^>/ {print $1}' genome.fa | sort | uniq -d finds them before they cause damage.

Assuming header text is preserved

samtools faidx and pyfaidx both use the first whitespace-delimited token as the record name. A header like >chr1 human chromosome 1 becomes chr1 in the index and in the BED. If you wanted the full description in the name column, you must extract it separately, and you must replace internal spaces with underscores before writing BED.

Irregular line wrapping

A FASTA whose sequence lines have different lengths within a record cannot be indexed directly. The error mentions a different line length in the named sequence. Re-wrap with seqkit seq -w 60 or awk before indexing.

Plain gzip instead of bgzip

samtools faidx works on uncompressed FASTA and on bgzip-compressed FASTA. A file produced by gzip will not index. Use bgzip genome.fa followed by samtools faidx genome.fa.gz. For large BED outputs, bgzip plus tabix -p bed gives you an indexed, randomly accessible file.

Chromosome naming mismatches

chr1 in the FASTA and 1 in an annotation file produce zero overlaps and a BED that looks structurally fine. Decide on one convention at the start and enforce it. This matters most when the BED will be intersected with a file from a different source.

Treating alignment coordinates as genome coordinates

An aligned FASTA is a matrix of alignment columns, not genome positions. A gap in one sequence consumes a column, so position 500 in the alignment is not position 500 in the genome. Converting an alignment to intervals requires mapping each ungapped segment back to reference coordinates. When gaps must be preserved, three-column BED is insufficient and BED12 blocks are the correct representation: a sequence ACGTACGTACGT aligned against ACGT--ACGT becomes chr1 0 12 seq2 0 + 0 12 0 2 4,4, 0,8,, which encodes two blocks of 4 bases separated by a 4-base gap. Reformatting sequence data is error-prone in exactly this way, which is why tools built for FASTA conversion put a review step in front of the output rather than writing blindly [6].

Confusing FASTA with FASTQ

They are not interchangeable and neither carries coordinates. A FASTQ record has a @ header, a sequence line, a + separator, and a quality line of the same length as the sequence. Converting FASTQ to FASTA throws away per-base quality. Converting FASTA to BED throws away the sequence entirely.

Limitations you cannot engineer away

BED cannot store sequence or quality. FASTA cannot store position. Any FASTA-to-BED conversion inherits whatever the headers claim, and if the headers are wrong the BED is wrong no matter how carefully you parse. A FASTA with no coordinates and no coordinate-bearing headers yields only whole-record intervals. That is the honest ceiling of the conversion, and it covers most of what people actually need.

Two more practical limits. Memory and speed: reading a 3 Gb FASTA with Biopython loads the sequences, while samtools faidx and pyfaidx read the index and seek, which is dramatically faster on repeat access. Parallelism: whole-record conversion is trivially parallel by chromosome, but header parsing is not, because you cannot split a FASTA by byte offset without knowing the record boundaries first.

Language choice also matters for teams. R users reach for Biostrings and rtracklayer, and the same file written by an R script and read by a Python script must round-trip without corruption. Cross-language read and write of a shared file format is an explicit design goal of modern bioinformatics packages precisely because mixed-language pipelines are normal [8]. Tab-delimited output is the safest interchange layer for that reason, since spreadsheet and database tools read it without special handling [9].

Frequently Asked Questions

Can I convert FASTA to BED with a single command?

Yes, for whole-chromosome intervals. samtools faidx genome.fa followed by awk 'BEGIN{OFS="\t"} {print $1, 0, $2}' genome.fa.fai > genome.bed is a complete two-command conversion. Any other scenario needs parsing logic.

Why is my BED end coordinate one less than the sequence length?

It should not be. A whole-record BED for a sequence of length L is [0, L), so the end column equals L exactly. If you see L - 1, you copied an inclusive end or subtracted one from the wrong column.

Does samtools faidx work on compressed FASTA?

It works on plain text and on bgzip-compressed FASTA with a bgzip-generated index. Files compressed with ordinary gzip cannot be indexed. Convert with bgzip first.

What is the difference between FASTA and FASTQ?

FASTQ stores a per-base quality score in a fourth line for every record. FASTA stores sequence only. Both are coordinate-free, so neither converts to BED without an external source of positions.

How do I make a BED from a FASTA with no coordinates in the headers?

You cannot make a positional BED from it. You can make a whole-record BED from the sequence lengths, which is the only coordinate information the file contains.

Why does bedtools reject my BED file?

The usual causes are a plain header row, space-separated columns, a non-integer score, a start coordinate greater than or equal to the end, or fewer than three columns. Run the NF and start-versus-end awk checks first.

Can BED represent a gapped alignment?

Only with BED12 block fields. Three-column BED describes one continuous interval. Use blockCount, blockSizes, and blockStarts to encode each ungapped segment as its own block.

Should the BED file be sorted?

Sort it if the receiving tool requires sorted input, which many interval operations do. sort -k1,1 -k2,2n produces chromosome-then-position order. Whole-chromosome BEDs are already sorted by construction if the FASTA records are in order.

Sources

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