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:
cd ~/bioinformatics_training/part7
ls -lh mapped.sorted.bam mapped.sorted.bam.baiIf 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.
samtools coverage mapped.sorted.bamThe 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.
samtools depth -aa mapped.sorted.bam > depth.tsv
head depth.tsv
wc -l depth.tsvdepth.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.
awk '$2>=1048 && $2<=1054' depth.tsvThis 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.
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.
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.tsvBoth 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:
awk '$3==0{n++} END{print n+0}' depth.tsv
awk '{sum+=$3} END{print sum/NR}' depth.tsvThe 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.
| Measure | Tested result | Meaning |
|---|---|---|
| Depth table | 2,400 rows | Every 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 depth | 11.6667× | An average across all 2,400 positions, including zeros. |
| Zero-depth positions | 350 bases: 1051–1400 | The 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:
sudo apt update
sudo apt install minimap2 samtools -y
minimap2 --version
samtools --versionCreate a new folder and open it in Windows Explorer:
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:
minimap2 -ax sr reference.fasta reads.fastq.gz > mapped.sam
samtools sort -o mapped.sorted.bam mapped.sam
samtools index mapped.sorted.bamContinue 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:
bash commands.shThe 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.
- What mean depth does
samtools coveragereport? - How many reference positions have zero depth?
- What percentage reaches at least 1×?
- What percentage reaches at least 10×?
- 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.