SAM, BAM & CRAM

After sequencing you have millions of short reads in FASTQ files. The next step is to find where each read came from in the genome, a process called alignment or mapping. An aligner such as BWA or STAR takes each read, finds the best matching position on a reference genome, and writes an alignment file. SAM and BAM are the standard formats for that result.

SAM format structure

SAM stands for Sequence Alignment/Map. It is a plain text format with two sections. Header lines start with @ and describe the reference genome, the aligner and other metadata. Alignment lines follow, one line per read and its mapped position.

@HD  VN:1.6  SO:coordinate
@SQ  SN:chr1  LN:248956422
@SQ  SN:chr2  LN:242193529
@PG  ID:bwa  PN:bwa  VN:0.7.17
SRR123456.1  99  chr1  10050  60  150M  =  10200  300  ATCGATCG...  IIIIIIII...
SRR123456.1  147  chr1  10200  60  150M  =  10050  -300  GCTAGCTA...  IIIIIIII...

The @HD line gives the format version and the sort order. Each @SQ line describes one chromosome or contig in the reference. The @PG line records which program created the file.

The 11 mandatory SAM fields

Every alignment line has exactly 11 tab-separated fields. Optional fields may follow.

Column	Name	Description
1	QNAME	Read name. Matches the FASTQ header.
2	FLAG	Bitwise flag encoding read properties.
3	RNAME	Reference sequence name (chromosome).
4	POS	1-based leftmost mapping position.
5	MAPQ	Mapping quality score.
6	CIGAR	How the read aligns to the reference.
7	RNEXT	Reference name of the mate read. = means same chromosome.
8	PNEXT	Position of the mate read.
9	TLEN	Template length (insert size).
10	SEQ	Read sequence.
11	QUAL	Base quality scores in Phred+33, same as FASTQ.

The FLAG field

The FLAG is a single number that encodes several yes/no properties as bits.

Flag value	Meaning
1	Read is paired
4	Read is unmapped
8	Mate is unmapped
16	Read is on the reverse strand
64	First read in the pair (R1)
128	Second read in the pair (R2)
256	Secondary alignment
1024	PCR or optical duplicate

The bits combine by addition. A FLAG of 99 decomposes as 1 + 2 + 32 + 64: the read is paired, mapped in a proper pair, its mate is on the reverse strand, and it is the first read of the pair. You never decode these by hand. samtools flags translates for you.

samtools flags 99
# 0x63  99  PAIRED,PROPER_PAIR,MREVERSE,READ1

The MAPQ field

MAPQ estimates the probability that the read is mapped to the wrong position, on the same Phred scale as base qualities. A MAPQ of 60 means the aligner is very confident. A MAPQ of 0 means the read maps equally well to several locations. A common filter keeps only reads with MAPQ of 30 or more.

The CIGAR string

The CIGAR string describes how the read aligns to the reference in a compact notation.

Operation	Meaning
M	Alignment match (a match or a mismatch)
I	Insertion to the reference
D	Deletion from the reference
N	Skipped region (an intron in RNA-seq)
S	Soft clipping (bases in SEQ that did not align)

150M means all 150 bases align contiguously. 75M1I74M means 75 bases match, then a one-base insertion, then 74 more matches. 50M1000N100M means 50 bases match, then a 1000-base intron, then 100 bases match: the signature of a spliced RNA-seq read.

BAM versus SAM versus CRAM

SAM is human-readable but enormous, and a whole-genome SAM can reach hundreds of gigabytes. BAM is the binary, compressed version of SAM. It holds exactly the same information in a third to a fifth of the space, and it is the working format for almost every tool. CRAM goes further with reference-based compression: it stores only the differences between each read and the reference genome, which makes files 30 to 50 percent smaller than BAM. The trade-off is that you need the reference genome at hand to decompress a CRAM.

Tip

Always work with BAM, not SAM. Convert SAM to BAM right after alignment, and delete the SAM file.

# Convert SAM to BAM
samtools view -bS alignment.sam > alignment.bam

# Convert BAM to CRAM for long-term storage
samtools view -C -T reference.fa alignment.bam > alignment.cram

Sorting and indexing

Aligners write reads in the order they appear in the FASTQ file. Most downstream tools need coordinate-sorted BAM, where reads are ordered by chromosome and position, and an index file (.bai) that maps genomic coordinates to byte offsets so a tool can jump straight to a region.

# Sort by coordinate, then index
samtools sort alignment.bam -o alignment.sorted.bam
samtools index alignment.sorted.bam
# Creates alignment.sorted.bam.bai

Caution

You must sort before you index. An unsorted BAM cannot be indexed, and tools fail silently or return wrong results on an unsorted file.

Some tools want name-sorted input instead. featureCounts and htseq-count can count paired-end reads from a name-sorted BAM. Check the documentation of each tool.

Key samtools commands

samtools is the toolkit for BAM files. These are the commands you will use most.

View alignments

# View the header
samtools view -H alignment.sorted.bam

# View alignments in a specific region
samtools view alignment.sorted.bam chr1:10000-20000

# View only mapped reads
samtools view -F 4 alignment.sorted.bam

# View only unmapped reads
samtools view -f 4 alignment.sorted.bam

The -f flag keeps reads that carry the given flag bits. The -F flag removes them.

Alignment statistics

samtools flagstat alignment.sorted.bam

50000000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
2500000 + 0 duplicates
48000000 + 0 mapped (96.00% : N/A)
50000000 + 0 paired in sequencing
25000000 + 0 read1
25000000 + 0 read2
46000000 + 0 properly paired (92.00% : N/A)

Check two numbers first. The mapping rate should exceed 80 percent when the reference matches the sample. The properly paired rate should sit close to the mapping rate; a wide gap points to insert size problems or chimeric fragments.

Per-chromosome counts and depth

# One line per chromosome: name, length, mapped reads, unmapped reads
samtools idxstats alignment.sorted.bam

# Average depth across the genome
samtools depth alignment.sorted.bam | awk '{sum+=$3} END {print sum/NR}'

# Depth in one region
samtools depth -r chr1:10000-20000 alignment.sorted.bam

What to check in a new BAM file

Mapping rate. Expect 90 percent or more for DNA-seq and 70 to 90 percent for RNA-seq against a well-matched reference. Below that, suspect contamination or the wrong reference.

Duplicate rate. PCR duplicates are reads that came from the same DNA fragment, an artefact of library preparation. A duplicate rate above 30 percent suggests low library complexity. Mark them with picard MarkDuplicates or samtools markdup.

Insert size distribution. For paired-end data the insert size should form a tight distribution around the expected fragment length. A broad distribution or several peaks indicate library preparation problems.

Mapping quality distribution. Most reads should have high MAPQ. A large block of MAPQ 0 reads means the sample maps poorly to the reference or carries a lot of repetitive sequence.

Summary

• SAM is the text form of read alignments, BAM its compressed binary form, and CRAM a smaller reference-based form.
• Each alignment line has 11 fields. FLAG, MAPQ and CIGAR carry the most information.
• Always work with sorted, indexed BAM files.
• Use samtools flagstat for the mapping and duplicate rates, and samtools view with flag filters to extract read subsets.
• Check mapping rate, duplicate rate, insert size and mapping quality on every new BAM file.