Where do our reads belong?

Step 01 · Read placement

A reference gives the reads a coordinate system.

In Part 6, we considered whether reads needed cleaning. Now we ask where reads align on a chosen reference sequence. A read may align on either strand, contain a mismatch, or have no suitable alignment.

The reference is the comparison target you select. Alignment to it does not by itself prove the sample’s organism identity. Likewise, a high-quality read can remain unmapped when it does not match the reference under the selected settings.

1. 1
   Choose a comparison

   Use the supplied synthetic reference and reads.

2. 2
   Place the reads

   Record alignment positions, strand and differences.

3. 3
   Review the evidence

   Inspect mapped records and consider what remains unknown.

Step 02 · Follow the information

From FASTA and FASTQ to SAM, BAM and BAI.

reference.fasta

The selected reference DNA sequence, used as the coordinate system.

reads.fastq.gz

Compressed reads with base-call quality scores.

mapped.sam

Text-based alignment records produced by minimap2.

mapped.sorted.bam

Compressed binary alignments sorted by reference coordinate.

mapped.sorted.bam.bai

An index for efficient access to positions within the BAM.

Sorting changes the order of alignment records; it does not change the DNA sequences. The BAI is an index, rather than a second set of reads.

Step 03 · Ubuntu terminal

Install minimap2 and SAMtools once.

Open the Ubuntu environment you set up in Part 1 and run:

sudo apt update
sudo apt install minimap2 samtools -y
minimap2 --version
samtools --version

Your Ubuntu password may be requested. If both tools are already installed, check their versions and continue.

Step 04 · Keep the files together

Prepare a separate Part 7 folder.

mkdir -p ~/bioinformatics_training/part7
cd ~/bioinformatics_training/part7
explorer.exe .

In the Windows Explorer window, copy reference.fasta and reads.fastq.gz from the extracted practice ZIP into this folder. Return to Ubuntu and check:

ls -lh reference.fasta reads.fastq.gz

The reference is 2,400 bases long. The FASTQ contains 300 reads of 100 bases each: 280 were constructed from two reference flanks and 20 are unrelated synthetic background reads. Some reference-derived reads are reverse-complemented or contain one substitution. These are data-design counts; the aligner’s output still needs checking.

Step 05 · Create the alignment

Map the single-end short reads.

minimap2 -ax sr reference.fasta reads.fastq.gz > mapped.sam

-x sr selects minimap2’s short-read preset; -a requests SAM output. Together they appear as -ax sr. The > redirects the alignment output into mapped.sam. Progress messages may still appear in the terminal.

This exercise uses single-end short reads. Different sequencing data can require a different preset. Running this command again replaces the existing mapped.sam in the current folder.

Step 06 · Organise the alignments

Sort the BAM and create its index.

samtools sort -o mapped.sorted.bam mapped.sam
samtools index mapped.sorted.bam

The first command reads the SAM and writes a coordinate-sorted BAM. The second creates mapped.sorted.bam.bai. Both commands can finish without printing a success message.

Keep the reference, reads, sorted BAM and BAI together: you will use these files again for the next exercise.

Step 07 · Check the result

How many alignment records are mapped?

samtools flagstat mapped.sorted.bam
samtools view -c -F 4 mapped.sorted.bam

flagstat summarises alignment records and their flags. In view -c -F 4, -c counts records and -F 4 excludes records flagged as unmapped. Secondary or supplementary alignments can make record counts differ from read counts in other datasets.

For this tested dataset there are no secondary or supplementary records, so the mapped record count also equals the number of mapped reads.

What happened in the tested exercise?

The supplied commands were run in Ubuntu 22.04 with minimap2 2.24-r1122 and SAMtools 1.13. All commands completed successfully, and the sorted BAM passed an integrity check.

On a small screen, scroll the table sideways.

Exact outputs can vary with software versions or settings. Compare your results with these observations and investigate differences.

Step 08 · Ten-minute challenge

Explain the result in your own words.

1. Run mapping, sorting and indexing from your Part 7 folder.
2. Check the mapping summary and mapped record count.
3. Explain what information is stored in FASTQ, SAM, BAM and BAI.
4. Why does 93.33% mapping not establish complete reference coverage?
5. Why does an unmapped read not automatically mean failed sequencing?

If you get stuck

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

Further reading: the official minimap2 documentation and the SAMtools guides for sort, index and view.