How to Convert BED to SAM: A Guide

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

By the end of this guide you will have a SAM file built from a BED interval file, with a header samtools accepts, 1-based positions that match the original intervals, and CIGAR strings you can defend in a code review. You will also know exactly which BED columns survive the trip, which ones are silently dropped, and when the conversion is the wrong idea in the first place.

Have these on hand before you start:

• A BED file with at least three columns, tab-separated, sorted by chromosome and start.
• The reference FASTA for the assembly the BED was made against, for example GRCh38.
• bedtools version 2.25 or newer and samtools version 1.10 or newer. Both install with conda install -c bioconda bedtools samtools.
• Optional: pysam if you want the Python route, and a POSIX awk if you want the zero-dependency route.

A quick word on expectations. BED describes a region. SAM describes a read that has been placed on a reference. Moving data from the first model to the second creates a synthetic alignment record. That record is useful for tools that only read SAM or BAM, and it is useless for anything that needs real bases or base qualities. Format conversion is a re-expression of one data model in another, and the gap between the two models is exactly where information disappears [1]. Keep that in mind and the rest of this guide will make sense.

What BED and SAM Actually Store

The BED side of the equation

BED is a plain text, tab-delimited format defined by the UCSC Genome Browser group. It has three required columns and up to nine optional ones. Column 2 is zero-based and column 3 is exclusive, which is the half-open convention.

A three-column file is valid BED. Most peak callers emit at least BED6. Splice-aware tools emit BED12. You need to know which one you have before you pick a conversion command, because the block columns change the CIGAR string completely.

The SAM side of the equation

SAM is the tab-delimited alignment format maintained alongside BAM and CRAM. Every alignment line has 11 mandatory fields, in this order.

The header matters as much as the alignment lines. A SAM file needs an @HD line recording the version and sort order, and one @SQ line per reference sequence with SN (the name) and LN (the length). BED carries no chromosome lengths, so those LN values have to come from somewhere else. That is the single most common reason a BED to SAM conversion fails downstream.

Why POS is one higher than chromStart

BED is zero-based. SAM is one-based. The first base of a chromosome is position 0 in BED and position 1 in SAM. So an interval that starts at BED column 2 value 1000 begins at SAM POS 1001. This is not a formatting quirk. It is a different counting convention, and getting it wrong shifts every interval in the file by one base.

SAM does not store an end coordinate. The end is implied by POS plus the number of reference bases the CIGAR string consumes, minus one. A 500M CIGAR at POS 1001 spans 1001 through 1500, which corresponds to a BED interval of 1000 to 1500. The arithmetic closes cleanly only if you add one in the right place.

What is lost, and why it cannot be recovered

Two losses are unavoidable. First, SEQ and QUAL become *. The BED file never held sequence or base qualities, so no conversion tool can invent them. Second, the record has no alignment evidence behind it. There is no edit distance, no MD tag, no insert size, no mate information. Any converter that claims to preserve those is fabricating data.

A conversion routine that is not tested against these gaps will fail quietly. Libraries that handle many formats treat comprehensive testing and round-trip validation as part of the deliverable rather than an afterthought [2].

Before You Run Anything

Confirm the BED is really BED

Files with a .bed extension are not always BED. They can be BED-like output from a tool with a different column order, or a space-delimited table that was renamed. Parsing a file whose grammar is not what you assumed is a distinct problem from converting it, and it needs to be solved first [3].

head -n 2 regions.bed | cat -A | head -n 2

cat -A shows ^I for tabs and $ for line ends. If you see runs of spaces instead of ^I, the file is not tab-delimited and every downstream step will misread the columns.

Check that the field count is consistent.

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

You want one or two distinct counts, typically 6 and 6, or 12 and 12. A count of 1 means the delimiter is wrong.

Strip UCSC track and browser lines, and comment lines.

grep -v -E '^(track|browser|#)' regions.bed > regions.clean.bed
wc -l regions.clean.bed

That wc -l number is your expected SAM record count. Write it down.

Decide on a name for every interval

SAM QNAME cannot be empty, cannot contain whitespace, and is limited to 254 characters. BED files from peak callers usually have usable names. BED files from coverage tools often do not.

Check for missing names.

awk -F'\t' 'NF < 4 || $4 == "" || $4 == "." {c++} END {print c+0, "records need a name"}' regions.clean.bed

If that number is large, plan to synthesize names like interval_1, interval_2. Do not leave QNAME as *. Many tools key on read name, and duplicates cause trouble later.

Pick BED6 or BED12 deliberately

If your file has 12 columns and every row has a valid blockCount that matches the number of entries in blockSizes and blockStarts, you want the BED12 route. If the block columns are present but empty or inconsistent, you want the BED6 route and you should accept that the introns will be filled in. Running a BED12 conversion on malformed block columns is a reliable way to produce a CIGAR string that samtools rejects.

How to Convert BED to SAM Step by Step

Step 1: Build the genome file that supplies @SQ

The genome file is a two-column, tab-delimited table of sequence name and length. There is no header. If you already have the reference FASTA, samtools can index it and you can cut the lengths straight out of the .fai.

samtools faidx GRCh38.fa
cut -f1,2 GRCh38.fa.fai > GRCh38.genome
head -n 3 GRCh38.genome

Expected output, for GRCh38:

chr1	248956422
chr2	242193529
chr3	198295559

The .fai file has five columns, and cut -f1,2 keeps the name and the length. If you do not have the FASTA, the UCSC hg38.chrom.sizes file is the same thing with the hg38 naming convention (chr1, and so on). Do not mix an Ensembl-style BED (1, 2, X) with a UCSC-style genome file (chr1, chr2, chrX). The names must match character for character.

Step 2: Sort the BED

Most converters require the BED to be grouped by chromosome, with records in increasing start order within each chromosome. Sorting lexicographically works because each chromosome stays contiguous. Sorting in the assembly's own chromosome order is tidier.

bedtools sort -i regions.clean.bed -g GRCh38.genome > regions.sorted.bed

If you do not want a bedtools dependency at this stage, plain sort is enough:

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

The -k2,2n numeric sort is essential. Without the n, start position 1000 sorts before start position 200, and the file is out of order in a way that only shows up on chromosome 2.

Step 3: Run bedtools bedtobam

There is no samtools bedtobam command. The conversion lives in bedtools, which writes BAM. You then expand the BAM to SAM text with samtools.

bedtools bedtobam \
  -i regions.sorted.bed \
  -g GRCh38.genome \
  -mapq 255 \
  > regions.bam

The three options that matter:

• -i is the input BED. It must be sorted.
• -g is the genome file. Without it, bedtools has no way to write @SQ lines and stops.
• -mapq sets the MAPQ value on every record. The default is 255, which is the SAM convention for "mapping quality unavailable". That is the honest value for a BED-derived record.

Add -bed12 when the input is BED12, and -ubam if you want uncompressed BAM, which is useful when you are about to stream the output straight into samtools.

Step 4: Expand the BAM to SAM

samtools view -h regions.bam > regions.sam

The -h flag includes the header. Without it you get alignment lines only, and a SAM file with no @SQ lines is rejected by almost everything.

Inspect the top of the file.

@HD	VN:1.6	SO:unknown
@SQ	SN:chr1	LN:248956422
@SQ	SN:chr2	LN:242193529
@SQ	SN:chr3	LN:198295559
peakA	0	chr1	1001	255	500M	*	0	0	*	*

The fields in a SAM file are tab-separated. They are shown with spaces above for readability. Your version string on the @HD line may differ, and the number of @SQ lines depends on how many sequences are in the genome file.

The first alignment line reads as follows: query name peakA, FLAG 0 (forward strand, mapped), reference chr1, position 1001, MAPQ 255, CIGAR 500M, no mate reference, mate position 0, template length 0, no sequence, no qualities.

Step 5: Work through a BED12 example

Take a two-exon transcript on the plus strand.

chr1	1000	1500	tx1	0	+	1000	1500	0	2	100,100	0,400

The blocks start at offsets 0 and 400 from chromStart, and each block is 100 bases. So exon one covers 1000 to 1100 and exon two covers 1400 to 1500. The intron between them is 300 bases.

bedtools bedtobam -i spliced.sorted.bed -g GRCh38.genome -bed12 -mapq 255 > spliced.bam
samtools view spliced.bam | head -n 1

tx1	0	chr1	1001	255	100M300N100M	*	0	0	*	*

The CIGAR reads as 100 matched bases, a 300-base reference skip, and 100 more matched bases. The reference span is 100 plus 300 plus 100, which is 500, matching chromEnd minus chromStart. The N operation is the standard way to encode an intron, and it is why splice-aware tools can read a BED12-derived SAM without special handling.

Note what the thickStart and thickEnd columns did. Nothing. They were dropped. A SAM record has no place to put coding region boundaries.

Step 6: Verify before you move on

samtools quickcheck -v regions.bam
samtools view -c regions.bam
wc -l < regions.sorted.bed

The two counts should match. If they do not, the converter skipped records, which usually means blank lines or a chromosome missing from the genome file.

Count the @SQ lines and compare with the genome file.

samtools view -H regions.sam | grep -c '^@SQ'
wc -l < GRCh38.genome

If your genome file has 25 sequences and your SAM has 24 @SQ lines, one sequence was dropped and every interval on it went with it.

Tips before you run the conversion on a new dataset

Run the pipeline on the first 100 lines first. Most coordinate and header problems surface immediately.

head -n 100 regions.sorted.bed > test.bed
bedtools bedtobam -i test.bed -g GRCh38.genome > test.bam
samtools view -h test.bam | head -n 5

Sort the test file also, because head does not preserve sort order in general even when the source is sorted. In practice it does for a sorted file, but check.

Set -mapq 0 instead of 255 if the synthetic records will pass through a variant caller or a coverage tool that filters on MAPQ. A MAPQ of 0 tells those tools to discard the record, which is often what you want, because the record has no alignment evidence behind it anyway.

Keep the BAM. A SAM text file is roughly five to ten times larger than the BAM it came from and takes longer to read. Store the BAM, and generate SAM only when a specific tool demands text.

Doing the Same Job in Python with pysam

The shell route is fine for BED6 and BED12. Use Python when you need conditional logic, custom QNAMEs, or a header with extra @RG or @CO lines. Converters for other formats commonly augment the output beyond what the source document contained, and a SAM header is a good place for that kind of provenance [4].

import re
import pysam

OPS = {"M": 0, "I": 1, "D": 2, "N": 3, "S": 4, "H": 5, "=": 7, "X": 8}

def cigar_tuples(cigar_string):
    return [(OPS[op], int(n)) for n, op in re.findall(r"(\d+)([MIDNSHP=X])", cigar_string)]

def build_cigar(start, end, block_sizes=None, block_starts=None):
    if not block_sizes:
        return [(0, end - start)]
    sizes = [int(x) for x in block_sizes.split(",") if x]
    starts = [int(x) for x in block_starts.split(",") if x]
    ops = []
    for i, size in enumerate(sizes):
        if i:
            gap = starts[i] - (starts[i - 1] + sizes[i - 1])
            if gap > 0:
                ops.append((3, gap))
        ops.append((0, size))
    return ops

genome = {}
with open("GRCh38.genome") as fh:
    for line in fh:
        chrom, length = line.rstrip("\n").split("\t")[:2]
        genome[chrom] = int(length)

header = {
    "HD": {"VN": "1.6", "SO": "unknown"},
    "SQ": [{"SN": c, "LN": n} for c, n in genome.items()],
}

with pysam.AlignmentFile("regions.sam", "w", header=header) as out:
    with open("regions.sorted.bed") as fh:
        for n, line in enumerate(fh, 1):
            f = line.rstrip("\n").split("\t")
            chrom, start, end = f, int(f[5]), int(f[1])
            name = f[2] if len(f) > 3 and f[2] not in ("", ".") else f"interval_{n}"
            strand = f[6] if len(f) > 5 else "+"
            sizes = f[7] if len(f) >= 12 else None
            starts = f[8] if len(f) >= 12 else None
            rec = pysam.AlignedSegment()
            rec.query_name = name
            rec.flag = 16 if strand == "-" else 0
            rec.reference_id = out.get_tid(chrom)
            rec.reference_start = start
            rec.mapping_quality = 255
            rec.cigartuples = build_cigar(start, end, sizes, starts)
            rec.query_sequence = None
            rec.query_qualities = None
            out.write(rec)

Two things to notice. reference_start is zero-based in pysam, which matches BED directly, and pysam adds the one when it writes POS. That is why the script does not add one explicitly. Second, query_sequence = None writes * to the SEQ field, and samtools reads it back without complaint.

Writing with mode w produces SAM text. Writing with mode wb produces BAM. The header dictionary is the same either way.

A Dependency-Free awk Route

If you have a shell but not bedtools, a short awk program does the job for BED6. It reads the genome file first so the @SQ lines appear before any alignment line.

awk -F'\t' -v OFS='\t' '
BEGIN {
  print "@HD", "VN:1.6", "SO:unknown"
}
FNR==NR {
  print "@SQ", "SN:" $1, "LN:" $2
  next
}
{
  name = ($4 != "" && $4 != ".") ? $4 : "interval_" FNR
  flag = ($6 == "-") ? 16 : 0
  span = $3 - $2
  print name, flag, $1, $2 + 1, 255, span "M", "*", 0, 0, "*", "*"
}' GRCh38.genome regions.sorted.bed > regions.sam

Verify it like any other SAM.

samtools view -b regions.sam > regions.check.bam
samtools view -c regions.check.bam

This route cannot do BED12 splicing, does no name deduplication, and will produce a broken CIGAR (0M or a negative length) if any interval has zero or negative width. It is a good quick check and a poor production tool.

Option Reference

How to Check the Result Is Right

Five checks catch nearly every real defect.

Record count. samtools view -c regions.bam must equal the line count of the sorted BED. A shortfall means records were dropped.

Header completeness. samtools view -H regions.sam | grep '^@SQ' must list every sequence that appears in the BED, with the same lengths as the FASTA used downstream. A mismatch here does not break the file itself. It breaks the next tool that compares the header to a reference.

Position arithmetic. Pull the first few records and recompute by hand.

samtools view regions.sam | head -n 3 | cut -f1,4,6

For each line, POS minus 1 must equal the BED chromStart, and the reference-consuming length of the CIGAR must equal chromEnd minus chromStart.

Round trip. Convert back to BED and compare the fields that are supposed to survive. A converter that supports only one direction cannot be validated this way, which is one argument for preferring tools with both directions built in [9].

bedtools bamtobed -i regions.bam | cut -f1-4,6 > roundtrip.bed
cut -f1-4,6 regions.sorted.bed > original.bed
sort original.bed > original.sorted.bed
sort roundtrip.bed > roundtrip.sorted.bed
diff original.sorted.bed roundtrip.sorted.bed && echo "round trip clean"

Any line in that diff is a field that did not survive. Common offenders are names containing spaces and strand values other than + or -.

Indexability. If downstream tools need an index, the records must be coordinate sorted and the header must say so.

samtools sort -o regions.sorted.bam regions.bam
samtools index regions.sorted.bam

The sort step rewrites the @HD line with SO:coordinate. That is safer than editing the header by hand and it costs seconds.

Common Mistakes and Limitations

Leaving POS equal to chromStart. This is the most common error and the hardest to notice. Every interval shifts one base to the left. A variant at the first base of a target region falls outside the window. Peak overlaps shrink by one base at the left edge. The fix is a single plus one, either in your script or by delegating to a tool that does it for you.

Writing a SAM file with no @SQ lines. samtools will refuse to read the alignments, and the error message points at the header rather than at your coordinates. Always include the genome file.

Using a genome file from the wrong assembly. The conversion succeeds. The file looks fine. Then a tool compares the LN value to the reference FASTA and stops. Worse, if the two assemblies differ in length only slightly, nothing stops and you get quietly wrong coordinates. Match the BED, the genome file, and the FASTA to the same assembly release.

Treating the output as a real alignment. A BED-derived SAM record has no bases, no base qualities, and no alignment score. Do not feed it to a variant caller, a read realigner, or anything that computes a likelihood from base qualities. Those tools either refuse the file or produce nonsense.

Zero-length intervals. BED permits chromStart equal to chromEnd for an insertion point. SAM does not have a way to express a zero-length alignment because a CIGAR string must consume at least one reference base. Widen these to one base before conversion, or keep them out of the BED entirely.

Malformed BED12 blocks. When blockCount does not match the number of entries in blockSizes, or when the last block extends past chromEnd, the CIGAR string comes out wrong. Validate the block columns before you run -bed12.

Silent loss of score and thick coordinates. No error is raised because nothing is broken. The columns simply do not exist in the target format. If you need score downstream, keep a side table keyed on QNAME, or keep the original BED.

Duplicate QNAMEs. Mark-duplicate tools and name-sorted operations assume read names identify records. Peak names in a BED are frequently duplicated across chromosomes. Deduplicate or suffix names before conversion.

Whitespace in names. QNAME cannot contain spaces. A BED name column of peak 1 produces a SAM line with an extra field, and every column after it is misaligned. Replace spaces with underscores before converting.

Assuming conversion is always the right answer. Automatic conversion between file formats can lose data, and adopting the format a tool already reads is often more practical than converting between formats [10]. IGV reads BED directly. bedtools genomecov -bg -i regions.bed -g genome.txt computes coverage from BED without touching SAM. Convert only when the downstream tool genuinely requires SAM or BAM.

Where BED, SAM, FASTA, and FASTQ Fit Together

These four formats answer different questions, and confusing them is the root of most BED to SAM mistakes.

FASTA holds named sequences with no quality information. A header line beginning with > is followed by the residues. It is the format of reference genomes and protein sequences.

FASTQ holds named sequences with one quality score per base. A header line beginning with @ is followed by the sequence, a + separator line, and an ASCII-encoded quality string of exactly the same length as the sequence. It is the format that sequencers produce.

SAM holds alignments of sequences to a reference. BED holds intervals on a reference.

The practical consequence is that BED to SAM and FASTQ to SAM are different operations. FASTQ to SAM is alignment, and it requires an aligner such as BWA or minimap2. BED to SAM is record construction, and it requires only a genome file. One produces data. The other produces a container.

Converting between sequence representations is a well-studied problem in its own right because the formats encode different levels of confidence about each base [11]. That is exactly the distinction that matters here. A FASTQ record carries per-base confidence. A BED interval carries none. When you write * into the QUAL field, you are stating that the confidence is unavailable, and you are correct to do so.

If you genuinely need SEQ filled in, bedtools getfasta can pull the reference sequence for each interval, and the -s flag reverse complements intervals on the minus strand so the sequence matches the FLAG you set.

bedtools getfasta -fi GRCh38.fa -bed regions.sorted.bed -s -name > regions.fa

You would then have to merge those sequences into the SAM records yourself. QUAL remains * no matter what you do, because the reference genome has no quality scores.

Choosing between SAM text and BAM binary follows the same logic as any other format decision. Binary containers are smaller and faster to read, and format revisions that reduce file size are a recurring theme across bioinformatics and engineering disciplines alike [8]. High-throughput sequencing pipelines adopted compact binary representations for exactly this reason [12]. The performance gap is measurable. One comparison of mesh export formats recorded a median load time of 293 ms for native PLY files against a median export time of 51,639 ms for OBJ files [7]. Different file types, same lesson: the container is not free.

The practical rule is to keep BAM as the working format, generate SAM text only for tools that demand it, and convert back to BAM as soon as you are done.

Frequently Asked Questions

Do I need a reference genome to convert BED to SAM?

Yes, you need the chromosome lengths. SAM requires one @SQ line per reference with a name and a length, and BED carries neither. A two-column genome file built from a FASTA index is enough. The full FASTA is only needed if you want to fill in the SEQ field.

Why is the POS in my SAM one base higher than the BED start?

Because BED counts from zero and SAM counts from one. The first base of a chromosome is position 0 in BED and position 1 in SAM. POS equals chromStart plus 1, and that offset is correct, not a bug.

Can converting BED to SAM recover read sequences or base qualities?

No. BED contains neither. A converted record carries * in the SEQ and QUAL fields, which is the SAM convention for unavailable data. You can fill SEQ from a reference FASTA, but QUAL stays unavailable.

How do I convert a BED12 file that contains spliced blocks?

Add the -bed12 flag to bedtools bedtobam. The converter turns each block into a match operation and each gap into an N operation, so a two-exon feature becomes something like 100M300N100M. Validate the block columns first, because malformed blocks produce broken CIGAR strings.

What MAPQ should a converted BED interval get?

255, which means mapping quality unavailable. That is the bedtools default and the honest value, because there is no alignment behind the record. Use 0 instead if you want downstream tools to filter the record out.

Can I convert SAM back to BED?

Yes. bedtools bamtobed -i file.bam returns chrom, start, end, name, score, and strand. Anything that was dropped on the way in stays dropped, so a round trip gives you back the six core columns, not the original twelve-column file.

Why does bedtools reject my BED file as unsorted?

Converters require records grouped by chromosome with increasing start positions. A BED sorted by string value rather than numerically puts start 1000 before start 200. Sort with sort -k1,1 -k2,2n or bedtools sort -g genome.txt and try again.

Should I keep the SAM file or the BAM file?

Keep the BAM. SAM text is larger and slower to read, and the two are interchangeable through samtools view. Generate SAM only when a tool specifically requires text, then convert back.

Sources

1. Development of an Algorithm for Converting JSON Formats to XML by Forming its File Data Structure
2. IOData: A python library for reading, writing, and converting computational chemistry file formats and generating input files.
3. Parsing without a grammar: making sense of unknown file formats
4. KEGGtranslator: visualizing and converting the KEGG PATHWAY database to various formats
5. Converting Between PLY and Ballistic Research Laboratory-Computer-Aided Design (BRL-CAD) File Formats
6. High Accuracy Converting between Different Coordinate Systems by Means of Drawing Interchange File(DXF)
7. Mesh quality comparison of common three file formats derived from intraoral scans
8. File formats for curvilinear multi-beam writing of 193i and EUV masks
9. cd2sbgnml: bidirectional conversion between CellDesigner and SBGN formats.
10. Qualitative study of the incompatibility of indoor map file formats with location software applications
11. Converting between Various Sequence Representations
12. CoCo: An application to store High-Throughput Sequencing data in compact text and binary file formats