SH•Syed Adnan HaiderClinical & pathogen genomicsGet in touch
Menu

FASTQ to Phylogeny · Part 8 · Free beginner exercise

How much of the genome is covered?

From mapped reads to evidence across a reference.

Use SAMtools to measure per-base depth and coverage breadth, find a region without reads, and see how the answer changes with your depth threshold.

Join the Part 8 discussion on LinkedIn ↗
The same 2,400-base referenceContinue with your Part 7 BAMCompare ≥1× and ≥10×Depth thresholds change breadthNo new tool after Part 7Use your existing SAMtools installation

Everything for Part 8

Download. Measure. Interpret.

Get the eleven-slide lesson and the practice ZIP. If you completed Part 7, you can continue directly with the BAM and index you already made.

Practice ZIP: about 9 KB · Slides: eleven pages, about 100 KB

Inside the ZIP: reference.fasta, reads.fastq.gz, README.txt, DESIGN_NOTES.txt and commands.sh.

The reference and reads are identical to Part 7 and entirely synthetic. Extract the outer ZIP; keep the FASTQ compressed. The README explains both the continuation route and a fresh start from the supplied inputs.

Step 01 · Three different questions

Mapping, depth and breadth.

Mapping rate
What fraction of reads aligned to the selected reference?
Read depth
How many aligned read bases cover one reference position?
Coverage breadth
What fraction of reference positions meet a specified minimum depth?

A dataset can have a high mapping rate while reads cluster in a few regions. Mean depth is an average across the reported region, including zeros; it can hide positions with little or no support.

Breadth (%) = positions meeting the depth threshold ÷ total reference positions × 100. Always state the threshold as well as the percentage.

Step 02 · Recommended route

Use the alignment from Part 7.

Open Ubuntu and return to the folder where you completed reference mapping:

Ubuntu terminal · Continue from Part 7
cd ~/bioinformatics_training/part7
ls -lh mapped.sorted.bam mapped.sorted.bam.bai

If both files are present, continue below. If you have not completed Part 7, use the fresh-start instructions with the practice ZIP.

Step 03 · Reference overview

Ask SAMtools for the summary.

Ubuntu terminal · Summarise coverage
samtools coverage mapped.sorted.bam

The output has a header and one row for this single reference. Focus on three columns:

covbases
Number of reference bases with depth of at least 1×.
coverage
Breadth at 1×, expressed as a percentage.
meandepth
Average read depth across the reported region.

In this exercise the reported region is the full 2,400-base reference. The summary does not show whether support is evenly distributed.

Step 04 · Keep zero-depth positions

Look at every reference position.

Ubuntu terminal · Inspect every reference position
samtools depth -aa mapped.sorted.bam > depth.tsv
head depth.tsv
wc -l depth.tsv

depth.tsv has three tab-separated fields: reference name, one-based position and read depth. It has no header in this command. For this one-reference dataset, wc -l should report 2,400 rows.

-aa includes positions with zero depth and unused reference sequences. Keeping those positions is essential when the entire reference is your denominator. A table containing only covered positions would overstate breadth.

Step 05 · Inspect the missing middle

Find where the evidence stops.

Ubuntu terminal · Inspect the start of the gap
awk '$2>=1048 && $2<=1054' depth.tsv

This prints positions 1048–1054. In the tested output, depth becomes zero at position 1051. The complete designed gap extends from 1051 through 1400: 350 reference positions without read support.

A gap is a clue, not a diagnosis.

For these teaching data, the gap was deliberately built in. In real data, sample and library effects, reference choice, alignment and filtering can all contribute. A coverage pattern alone does not identify the cause.

Step 06 · State the threshold

How much reaches 1× or 10×?

The first command calculates breadth at ≥1×; the second calculates breadth at ≥10×. Both use all rows in the depth table as the denominator.

Ubuntu terminal · Calculate breadth at two thresholds
awk '{n++; if($3>0)c++} END{print 100*c/n}' depth.tsv
awk '{n++; if($3>=10)c++} END{print 100*c/n}' depth.tsv

Both outputs are percentages of reference positions, rather than counts of reads. Requiring more depth produces equal or smaller breadth.

You can also count zero-depth positions and calculate mean depth from the same table:

Ubuntu terminal · Count zero-depth bases and calculate mean depth
awk '$3==0{n++} END{print n+0}' depth.tsv
awk '{sum+=$3} END{print sum/NR}' depth.tsv

The first output is a number of bases; the second is average depth across the entire reference.

What happened in the tested exercise?

The supplied workflow was tested in Ubuntu 22.04 with minimap2 2.24-r1122 and SAMtools 1.13. Both the step-by-step commands and optional script worked from fresh inputs and an existing Part 7 BAM.

Observed measurements for the supplied synthetic dataset.
MeasureTested resultMeaning
Depth table2,400 rowsEvery position in the reference is included.
Breadth at ≥1×2,050 / 2,400 = 85.4167%These positions have at least one aligned read base.
Breadth at ≥10×1,599 / 2,400 = 66.625%Fewer positions meet the higher threshold.
Mean depth11.6667×An average across all 2,400 positions, including zeros.
Zero-depth positions350 bases: 1051–1400The deliberately uncovered central interval.

On a small screen, scroll the table sideways.

SAMtools and awk may display different numbers of decimal places. Exact measurements can also change with tool versions, input alignments or filtering settings.

The 1× and 10× thresholds illustrate how breadth changes. Meeting either threshold alone does not establish the accuracy of a base call or suitability for a particular analysis.

Alternative route · From the practice ZIP

Starting without the Part 7 BAM?

Install the tools and recreate the alignment

Run these commands inside Ubuntu if the tools are not already installed:

Ubuntu terminal · Install and check tools
sudo apt update
sudo apt install minimap2 samtools -y
minimap2 --version
samtools --version

Create a new folder and open it in Windows Explorer:

Ubuntu terminal · Create a fresh Part 8 folder
mkdir -p ~/bioinformatics_training/part8
cd ~/bioinformatics_training/part8
explorer.exe .

Copy all five files from the extracted practice ZIP into this folder. Return to Ubuntu and recreate the alignment:

Ubuntu terminal · Recreate the Part 7 alignment
minimap2 -ax sr reference.fasta reads.fastq.gz > mapped.sam
samtools sort -o mapped.sorted.bam mapped.sam
samtools index mapped.sorted.bam

Continue with the coverage summary while staying in this Part 8 folder. Do not switch to the Part 7 folder for this route.

Use the optional commands.sh script

With commands.sh in your working folder, run:

Ubuntu terminal · Run the optional complete exercise
bash commands.sh

The script uses an existing mapped.sorted.bam if available; otherwise it maps the supplied reference and reads first. It writes coverage_summary.tsv and depth.tsv, then prints the breadth, mean depth and gap checks. Running it again replaces those two summary files.

Step 07 · Mini exercise

Interpret your coverage report.

  1. What mean depth does samtools coverage report?
  2. How many reference positions have zero depth?
  3. What percentage reaches at least 1×?
  4. What percentage reaches at least 10×?
  5. Why can a high mapping rate coexist with a 350-base gap?
Check your reasoning

The tested mean depth is about 11.67×, but it is not the depth at every position. The 350-base gap has no read evidence. Breadth falls from about 85.42% at 1× to 66.625% at 10×, showing why a mapping percentage or an average cannot describe the whole reference.

If you get stuck

Share the command, error or result you are unsure about in the Part 8 LinkedIn discussion.

My BAM file cannot be found.

Check pwd and ls -lh. Continue in the Part 7 folder if using your previous alignment, or stay in the Part 8 folder if you followed the fresh-start route.

depth.tsv has fewer than 2,400 rows.

Confirm that you used samtools depth -aa and the supplied 2,400-base reference alignment. Omitting zero-depth positions changes the denominator for breadth.

SAMtools does not recognise coverage.

Check samtools --version. This lesson was tested with SAMtools 1.13; the Ubuntu installation commands in the fresh-start instructions provide the required tools.

The output differs from the example.

Confirm the input files, reference length, command options and tool versions. Depth reflects the alignments and filters used. Keep the whole reference, including zero-depth positions, in the denominator for these calculations.

Further reading: the official SAMtools documentation for coverage and depth explains the output fields and filtering options.

Coming next · Part 9

From mapped reads to a consensus genome.

Which positions can we trust? Connect the alignment and coverage evidence to the next analysis step.