SAMtools, and the SAM/BAM/CRAM trio

One thing, three ways

SAM, BAM and CRAMare the same data

An aligner tells you where each read landed. That answer is written in one of three encodings, and SAMtools reads all of them.

SAM — plain text

Readable, greppable, and enormous. Useful for looking at, and almost never for keeping.

BAM — binary

The same records, compressed and indexable. The working format for essentially everything.

CRAM — reference-based

Stores only differences from a reference genome. Much smaller, but you must keep that reference.

The commands

Six subcommandsdo nearly everything

SAMtools has dozens. In practice a pipeline uses the same handful, in roughly the same order, every time.

align, sort, index

$ samtools sort -@ 8 -o aln.sorted.bam aln.bam

$ samtools index aln.sorted.bam

$ samtools flagstat aln.sorted.bam

viewConvert between the three formats, filter by flag or region, and count. The one you use most.

sortPut records in coordinate order. Required before indexing, and before most downstream tools.

indexWrite the .bai or .crai beside the file, so a region can be fetched without a full scan.

flagstat and statsHow many reads mapped, paired, duplicated. The first thing to run when a result looks odd.

merge and markdupCombine lanes into one file, then flag PCR duplicates so callers can ignore them.

mpileup and depthPer-base coverage, and the pile-up format that variant callers read as their input.

The rule

Sort first,then index

Almost every SAMtools problem is one of these two steps missing, and the error message rarely says so plainly.

A BAM straight out of an aligner is in read order, which is the order the sequencer produced them in. Nearly everything downstream wants coordinate order instead — sorted by chromosome, then position. Once sorted, an index lets a tool jump straight to a region without reading the file, which is what makes samtools view aln.bam chr1:1-1000 instant rather than a full scan.

Read orderOut of the alignerthe order the sequencer made them

By positionAfter sortchromosome, then coordinate

Random accessAfter indexa region without a full scan

IGV, variant callers and coverage tools all assume both steps have happened. If a genome browser shows you nothing, or a caller reports no reads in an obviously covered region, check for the .bai file before checking anything else.

Flags

The numberin the second column

SAM flags are a bitfield, which is why filtering reads feels cryptic until you see how it works.

It is a bitfield

One integer holding a dozen yes/no answers about how that read aligned.

-f keeps, -F drops

Filter on the bits rather than the number. -F 0x904 is the standard cleanup.

Do not decode by hand

The explain-flags page exists for a reason. Nobody memorises what 163 means.

Each bit means one thing: paired, properly paired, unmapped, reverse strand, secondary, duplicate. A flag of 99 is several of those at once. -f keeps reads with a bit set, -F excludes them, so -F 0x904 is the usual way to drop unmapped, secondary and supplementary alignments in one go.

On a machine

SAMtools isalready built

Nothing here needs compiling. Pick a machine, and these two commands put it on there.

hub://samtools

$ dxflow workflow create --identity samtools hub://samtools

$ dxflow workflow start samtools

Start a machine What SAMtools needs

Read it, then run itPrepaid and by the hour, on a machine that is yours about a minute after you ask.

Start a machine All notes