Bioinformatics File Formats Explained: FASTA, FASTQ, SAM/BAM, VCF, GTF and BED

Chapter 3 of 5: Alignment formats – SAM, BAM and CRAM

SAM, BAM and CRAM: alignment file formats

After alignment, each read gets a position on the reference. SAM (Sequence Alignment/Map) stores this as text. BAM is the same information in compressed binary form, and CRAM compresses further by storing only the differences from the reference (Li et al., 2009; SAM specification). Learn SAM, and you can read all three.

Step 3: Look at the header

A SAM file has a header (lines starting with @) and an alignment section (one line per alignment). Print just the header with -H:

samtools view -H toy.sam

The @SQ lines list every reference sequence with its name (SN) and length (LN). They match toy.fa.fai exactly, which is what you want. The @PG line was added by samtools itself: tools record each program run in the header, so a BAM carries its own history.

Step 4: The 11 mandatory columns

Every alignment line has at least 11 tab-separated columns, always in this order. Optional tags follow from column 12.

samtools view toy.sam | cut -f 1-6 | head -n 6

Note column 10: for reads that aligned to the reverse strand, SAM stores the reverse complement of the read, so that SEQ always matches the forward strand of the reference. The FLAG tells you which reads were flipped.

Step 5: Decode the FLAG

FLAG is a sum of powers of two. Each bit answers one yes/no question. Let samtools decode it:

samtools flags 163
samtools flags 83
samtools flags 4

So 163 = 128 + 32 + 2 + 1: paired, properly paired, mate on the reverse strand, and this is R2. Its mate, 83 = 64 + 16 + 2 + 1, is R1 on the reverse strand. Use -f to keep reads that have bits set, and -F to drop reads that have them:

samtools view -c toy.sam
samtools view -c -f 16 toy.sam
samtools view -c -F 0x904 toy.sam

• -c: count instead of printing.
• -f 16: only reads on the reverse strand.
• -F 0x904: drop unmapped (0x4), secondary (0x100) and supplementary (0x800) records, leaving one primary, mapped line per read. This is the usual way to count mapped reads.

# Is a read on the reverse strand? (pseudo-code)
is_reverse = (flag & 16) != 0
is_primary = (flag & 0x900) == 0      # neither secondary nor supplementary
is_mapped  = (flag & 4) == 0

# awk without bitwise operators: bit b is set if int(flag / b) % 2 == 1
awk '!/^@/ && int($2/16)%2==1' toy.sam    # reverse-strand records

Step 6: Read a CIGAR string

CIGAR is a list of length + operation pairs. Read r001 has 8M4I4M1D3M: 8 aligned bases, 4 inserted bases, 4 aligned, 1 deleted, 3 aligned.

Two sums check any alignment. Read length = sum of M I S = X = 8 + 4 + 4 + 3 = 19, the length of SEQ. Reference span = sum of M D N = X = 8 + 4 + 1 + 3 = 16, so r001 covers reference bases 7 to 22. The spliced read r004 (6M14N1I5M) jumps over 14 reference bases with N, exactly as RNA-seq reads jump over introns.

# Reference end of an alignment (pseudo-code)
end = POS - 1
for (length, op) in parse(CIGAR):
    if op in 'MDN=X':
        end += length
# end is 1-based and inclusive: r001 -> 7 - 1 + 16 = 22

Step 7: Convert to BAM, sort and index

Almost every downstream tool wants a coordinate-sorted, indexed BAM. One command sorts and writes BAM, another builds the index:

samtools sort -o toy.bam toy.sam
samtools index toy.bam
ls -l toy.sam toy.bam toy.bam.bai
samtools view toy.bam ref:30-40 | cut -f 1-6

• samtools sort -o toy.bam: sorts by reference and position and writes BAM (the format comes from the .bam name).
• samtools index: writes toy.bam.bai. The BAM must be coordinate-sorted first.
• samtools view toy.bam ref:30-40: prints only alignments that overlap bases 30–40 of ref (1-based region). This needs the index.

Notice r004: it starts at 16, but its 14N carries it past position 30, so it overlaps the region. Region queries use the full alignment span, not just POS.

samtools flagstat toy.bam
samtools idxstats toy.bam

idxstats reads the index and prints, per sequence: name, length, mapped reads and unmapped reads. On a real BAM this gives reads per chromosome in a second, which is a quick sex check (chrY) or contamination check (unexpected contigs).

Step 8: CRAM, the smallest alignment format

CRAM stores each read as differences from the reference, so it needs the same reference FASTA to write and to read.

samtools view -C -T toy.fa -o toy.cram toy.bam
samtools index toy.cram
samtools view -T toy.fa toy.cram | head -n 1

• -C: write CRAM.
• -T toy.fa: the reference used for compression and decompression.

Our toy CRAM (1,512 bytes) is actually bigger than the BAM (547 bytes), because CRAM has a fixed container overhead. Compression only pays off at real scale. On a simulated set of 100,000 read pairs (150 bp) against a 200 kb genome, we measured:

(Your numbers depend on your data.) Real sequencer data, with binned qualities and many identical bases, usually compresses better than our simulation, which used noisy, unbinned qualities.