The samtools Survival Kit

Before you start

You've done Read Alignment with BWA-MEM, so you have a BAM file to work on.
Comfortable running a few command-line lines. A free Google account to run the notebook.
Any term new? The Glossary has it.

Learning objectives

By the end of this lesson you will be able to: use the core samtools commands to look inside a BAM, sort reads into genome order, index the sorted BAM, summarize and pull a region, and check how deep the coverage is at a position.

Why a raw BAM isn't ready yet

When BWA finishes, the reads in your BAM file are in the order they came off the sequencer, which is essentially random. To analyze a genome you need to ask questions like "show me every read covering position 1,500 on chr1." Answering that on an unsorted file means scanning the whole thing every time. The fix is two quick steps: sort the reads into genome order, then build an index so any region can be found instantly. Almost every genomics tool expects a sorted, indexed BAM, so this is the universal next step after alignment.

Decode the jargonsamtools

samtools is the standard toolkit for working with SAM, BAM, and CRAM alignment files. It does conversions (SAM to BAM), sorting, indexing, filtering, statistics, and more. You call it as samtools <subcommand>, like samtools sort or samtools view.

▶ Run it in Google Colab

Colab is a free Linux machine in your browser, so samtools runs with nothing to install on your computer.

Open the notebook in Colab →

The five commands you'll actually use

Command	What it does
`samtools view`	Look inside a BAM, convert SAM to BAM, or filter reads.
`samtools sort`	Reorder reads by their genome position.
`samtools index`	Build the companion `.bai` index for fast region lookups.
`samtools flagstat`	One-glance summary: total, mapped, duplicate, paired reads.
`samtools depth`	Report how many reads cover each position (coverage).

1Look inside the BAM

BAM is binary, so you cannot just cat it. samtools view turns it back into readable text. Add -H to see just the header.

!samtools view -H aligned.bam        # the header
!samtools view aligned.bam | head -3   # first 3 alignments

2Sort into genome order

This is the key step. samtools sort reorders the reads so they run from the start of each chromosome to the end. We write the result to a new file.

!samtools sort aligned.bam -o aligned.sorted.bam

Decode the jargonCoordinate-sorted

A coordinate-sorted BAM has its reads ordered by chromosome and then by position along that chromosome. It's the order that lets tools (and the index) jump straight to any region instead of reading the whole file. When a tool complains it "expects a sorted BAM," this is what it means.

3Index the sorted BAM

Indexing creates a small companion file (aligned.sorted.bam.bai) that acts like a table of contents. With it, a tool can leap directly to "chr1:1500" without scanning everything before it. You must sort first; you cannot index an unsorted BAM.

!samtools index aligned.sorted.bam
!ls -lh aligned.sorted.bam*

4Summarize and pull a region

With a sorted, indexed BAM you get two superpowers. flagstat gives you the health summary, and view with a region string instantly returns just the reads over that spot.

!samtools flagstat aligned.sorted.bam
# Pull only the reads covering positions 1000-1100 on chr1:
!samtools view aligned.sorted.bam chr1:1000-1100 | wc -l

That region query is the payoff of indexing. On a real, multi-gigabyte BAM it returns in a fraction of a second instead of minutes.

⚠️ Sort before you index

Indexing only works on a coordinate-sorted BAM. If you try to index the raw output of the aligner, samtools will refuse. The reliable habit is always the same order: align, sort, index.

5How deep is the coverage?

Coverage (or depth) is how many reads overlap a given position. It's one of the most important quality numbers in genomics: variant callers trust a position more when many reads agree on it.

!samtools depth aligned.sorted.bam | head
# average depth across all covered positions:
!samtools depth aligned.sorted.bam | awk '{sum+=$3} END {print "mean depth:", sum/NR}'

Decode the jargonCoverage / depth

Coverage (often written like "30x") is the average number of reads stacked over each base. Higher coverage means more evidence at every position, so you can call variants with more confidence. Too little coverage and you can't tell a real variant from a sequencing error.

Where this leads

A sorted, indexed BAM with known coverage is exactly what a variant caller wants. In the next and final Track 2 lesson, you'll feed this file to a variant caller and produce a VCF, the list of every place your sample differs from the reference. That's the capstone that ties the whole workflow together.

🔶 Level up: mark duplicates

Real pipelines add one more step before variant calling: marking PCR duplicates (identical reads created during library prep, which can bias variant calls). Tools like samtools markdup or Picard's MarkDuplicates flag them so callers can ignore them. The principle is the same survival-kit thinking: clean the BAM before you trust it.

Check your understanding

Why must a BAM be sorted before it can be indexed?

The index is a table of contents that points to where each genome region sits in the file. That only makes sense if reads are stored in genome order. An unsorted BAM has reads scattered randomly, so there's nothing orderly for the index to point into, and samtools will refuse to index it.

What does a region query like chr1:1000-1100 rely on?

It relies on the BAM being coordinate-sorted and indexed. The index lets samtools jump straight to that slice of the file instead of scanning every read, which is what makes the query nearly instant on huge files.

You see a mean depth of 2x over your region. Is that enough to call variants confidently?

Usually not. At 2x, a single sequencing error can masquerade as a variant because so few reads vouch for each position. Most variant-calling work aims for much higher depth (often 30x or more for human germline work) so multiple reads independently confirm each call.

Which samtools subcommand reports how many reads cover each position (the coverage)?

samtools depth reports the number of reads overlapping each position. samtools sort reorders reads and samtools view looks inside or converts the file.

Why do real pipelines mark PCR duplicates before variant calling?

PCR and optical duplicates are copies of the same original fragment. Left in, they inflate the apparent evidence at a position and can cause false variant calls, so tools like samtools markdup or Picard flag them first.

Next in Track 2 · The capstone