An extra exercise for when you have finished Day 4: we will map the reads of sample 2 against the genome they come from, and measure how well they map and how deep the sequencing is.

General information:

• Raw reads: /shared/team/2026_training/week2/reads/lab/
• Your working directory: /shared/team/2026_training/week2/{firstname}
• As on Day 4, work only inside your own directory, and keep adding the commands and the answers to your notebook.md (create a new ## section for this exercise)
• Do not copy and paste blindly: read the output of every command

0. Start from your directory

Open a terminal and define the variables again (they are lost every time you open a new terminal):

cd /shared/team/2026_training/week2/{firstname}
export WORK=$PWD
export READS=/shared/team/2026_training/week2/reads
mkdir -p $WORK/mapping

All the output of this page will be saved in $WORK/mapping/.

1. A conda environment: mapping

Create a new environment called mapping with these tools:

mamba create -n mapping -c conda-forge -c bioconda bwa samtools bamtocov minimap2 ncbi-datasets-cli
conda activate mapping

Check that the tools work (note that bwa has no --version: it prints its version when run without arguments):

samtools --version | head -n 1
bamtocov --version
minimap2 --version
datasets --version
bwa 2>&1 | grep Version

Notebook: note the command used to create the environment.

2. Downloading the reference genome

Note that you can simply download the genome from NCBI and upload it to CLIMB. We see here how to do the same from the Command Line Interface

Sample 2 is Lambda DNA (the genome of the bacteriophage Escherichia phage Lambda, a classic control in sequencing). We will download its reference genome from NCBI using datasets.

Let’s search for it first, asking only for the NCBI RefSeq (the curated) assembly:

datasets summary genome taxon "Escherichia phage Lambda" \
   --assembly-source refseq --as-json-lines | head -c 600

The output is in JSON format, and it reports the accession of the assembly: GCF_000840245.1.

The accession of a genome is like its ID card: GCF_ is a RefSeq assembly, GCA_ is a GenBank assembly.

Create a directory for the genome, and download it:

mkdir -p $WORK/mapping/genome
cd $WORK/mapping/genome

datasets download genome accession GCF_000840245.1 --include genome --filename lambda.zip
unzip lambda.zip

The downloaded package contains several files: find the FASTA file of the genome (.fna):

find . -name "*.fna"

Copy it to a file with a simple name, lambda.fa, and have a look at it:

cp ncbi_dataset/data/GCF_000840245.1/*.fna lambda.fa
head -c 300 lambda.fa
grep ">" lambda.fa
samtools faidx lambda.fa
cat lambda.fa.fai

How many sequences are there? How long is the genome? (samtools faidx creates a small index, lambda.fa.fai, with one line per sequence: name, length…)

Notebook: note the accession, the name and the length of the genome.

3. Indexing the genome: bwa index

Before mapping, bwa needs to build an index of the reference (a data structure that allows searching it quickly). This has to be done only once for each genome:

cd $WORK/mapping/genome
bwa index lambda.fa
ls -lh

You will see new files next to lambda.fa, with extensions .amb, .ann, .bwt, .pac and .sa: together with lambda.fa they form the index, and they must stay in the same directory. The Lambda genome is small (about 48 Kbp) so it takes less than a second; for a human genome it takes about an hour.

4. Mapping sample 2: bwa mem

Sample 2 was sequenced in two pools, DTP1-2 and DTP2-2. We will start with the first one. Find its two files:

ls $READS/lab/PID-2959-DTP1-2_S*_R?_001.fastq.gz

The basic command for paired-end reads is:

bwa mem  reference.fa  reads_R1.fastq.gz  reads_R2.fastq.gz  > output.sam

The output is in SAM format: a text file with one line per read. It is very big, and slow to read, so we will never save it: instead we send it directly to samtools, using a pipe (|).

Making a sorted BAM

The BAM file is the compressed (binary) version of a SAM file, and it is the format used to store alignments. Most tools (including the coverage tools) also require it to be sorted by position on the reference and indexed.

samtools sort does the sorting and, writing the output to a .bam file, also the compression:

mkdir -p $WORK/mapping/bam

SAMPLE=DTP1-2

# Complete this command with the appropriate FILE PATHS
bwa mem -t 4 \
  $WORK/mapping/genome/lambda.fa \
  -Sample_R2_fastq- \
  -Sample_R2_fastq- \
  | samtools sort -@ 4 -o $WORK/mapping/bam/OUTPUTFILE.sorted.bam -

How it works:

• -t 4 and -@ 4 use 4 threads for bwa and samtools
• the \ at the end of a line means “the command continues on the next line”
• $(ls ...) finds the files, whose names contain a number we do not know in advance (_S7_)
• the pipe | sends the SAM output of bwa to the input of samtools sort, and the final - means “read from the pipe”
• this takes about 1 minute, and bwa will print a lot of messages: they are normal

Then index the BAM file:

samtools index $WORK/mapping/bam/${SAMPLE}.sorted.bam
ls -lh $WORK/mapping/bam/

You should see the BAM file (about 130 MB) and its index, DTP1-2.sorted.bam.bai. Each bam needs its own .bai next to it, with the same name.

Have a look inside the BAM file (samtools view converts it to text, and head shows the first lines):

samtools view -h $WORK/mapping/bam/${SAMPLE}.sorted.bam | head -n 10 | cut -c 1-120

The lines starting with @ are the header (the sequences of the reference, the program used), the other lines are the reads: read name, flag, reference, position, mapping quality, CIGAR…

Notebook: note the full command used, and the size of the BAM file.

5. How many reads are not mapped?

The quick summary: samtools flagstat

mkdir -p $WORK/mapping/stats
samtools flagstat $WORK/mapping/bam/${SAMPLE}.sorted.bam | tee $WORK/mapping/stats/${SAMPLE}.flagstat.txt

(tee prints the output on the screen and saves it in the file)

The first line is the total number of records, and the mapped line gives the number of reads aligned to the reference and, in brackets, the percentage.

Counting with flags

Each read in a BAM file has a flag: a number encoding its properties (paired, mapped, first in pair…). samtools view -c counts the reads, and we can select them by flag:

# all the reads
samtools view -c $WORK/mapping/bam/${SAMPLE}.sorted.bam

# unmapped reads
samtools view -c -f 4 $WORK/mapping/bam/${SAMPLE}.sorted.bam

# mapped reads
samtools view -c -F 4 $WORK/mapping/bam/${SAMPLE}.sorted.bam

The first number should be equal to the sum of the other two. Add the percentage of unmapped reads, using awk:

TOTAL=$(samtools view -c $WORK/mapping/bam/${SAMPLE}.sorted.bam)
UNMAPPED=$(samtools view -c -f 4 $WORK/mapping/bam/${SAMPLE}.sorted.bam)
echo "$UNMAPPED $TOTAL" | awk '{printf "%d unmapped reads out of %d (%.2f%%)\n", $1, $2, 100 * $1 / $2}'

You can also have a per-reference summary (reference name, length, mapped reads, unmapped reads with a mate mapped here), with:

samtools idxstats $WORK/mapping/bam/${SAMPLE}.sorted.bam

The last line (*) are the reads with no position at all.

You should find that about 99.7% of the reads map, and only ~11,000 reads (out of ~3.3 million) do not. Why do you think some reads do not map to the Lambda genome? (the sample was prepared in a laboratory, and sequenced on a machine that sequences many samples…)

Notebook: note the number and the percentage of unmapped reads.

Looking at the unmapped reads (optional)

Save the unmapped reads in FASTQ format, to inspect them (or to run Kraken2 on them!). Here we use the flag 12 (= 4 + 8), that selects the reads that are unmapped and have an unmapped mate, so that we obtain proper pairs (most of the other unmapped reads have a mate that maps, and they are a different story):

samtools fastq -f 12 -n \
  -1 $WORK/mapping/stats/${SAMPLE}.unmapped_R1.fastq \
  -2 $WORK/mapping/stats/${SAMPLE}.unmapped_R2.fastq \
  -0 /dev/null -s /dev/null \
  $WORK/mapping/bam/${SAMPLE}.sorted.bam

wc -l $WORK/mapping/stats/${SAMPLE}.unmapped_R?.fastq     # 4 lines per read

6. How deep is the coverage?

The coverage (or depth) of a position is the number of reads covering it. We are interested in:

• the breadth: what percentage of the genome is covered by at least one read
• the mean depth: how many times, on average, a base is sequenced

Quick summary: samtools coverage

mkdir -p $WORK/mapping/coverage
samtools coverage $WORK/mapping/bam/${SAMPLE}.sorted.bam | tee $WORK/mapping/coverage/${SAMPLE}.samtools-coverage.tsv | column -t -s $'\t'

For every reference sequence it reports: numreads, covbases (covered bases), coverage (percentage), meandepth, mean base quality and mean mapping quality.

There is also a text histogram, to see how the coverage is distributed along the genome:

samtools coverage -m $WORK/mapping/bam/${SAMPLE}.sorted.bam

What is the mean depth? Is the coverage uniform along the genome, or are there regions with less coverage?

Coverage per base: bamtocov

To see the depth at every position, we use bamtocov, which prints the regions with a constant coverage in BED format (reference, start, end, depth):

bamtocov $WORK/mapping/bam/${SAMPLE}.sorted.bam > $WORK/mapping/coverage/${SAMPLE}.bamtocov.bed
head $WORK/mapping/coverage/${SAMPLE}.bamtocov.bed
wc -l $WORK/mapping/coverage/${SAMPLE}.bamtocov.bed

Each line is a stretch of the genome with a given depth. We can find the minimum and the maximum depth, for example, using sort:

# Highest depth
sort -k4,4nr $WORK/mapping/coverage/${SAMPLE}.bamtocov.bed | head -n 1

# Lowest depth (ignoring the first line, that is the position 0)
tail -n +2 $WORK/mapping/coverage/${SAMPLE}.bamtocov.bed | sort -k4,4n | head -n 1

If you do not need the per-base table, bamtocov can also produce only a short report:

bamtocov --skip-output --report $WORK/mapping/coverage/${SAMPLE}.bamtocov.report.txt $WORK/mapping/bam/${SAMPLE}.sorted.bam
cat $WORK/mapping/coverage/${SAMPLE}.bamtocov.report.txt

Notebook: note the commands, and answer:

1. What percentage of the Lambda genome is covered?
2. What is the mean depth? Is it a very high or a very low coverage for a 48 Kbp genome?
3. Which are the regions with the lowest and the highest depth? (you can search them with sort, or open the .bed file in VisiData)

7. Do it again for the second pool

Repeat the mapping, the flagstat and the coverage for the pool DTP2-2: change the SAMPLE variable and re-run the commands (use the up arrow, do not retype everything!).

SAMPLE=DTP2-2

Does the second pool have a similar percentage of mapped reads and a similar mean depth?

At the end you should have this:

{firstname}/
└── mapping/
    ├── genome/
    │   ├── lambda.fa
    │   └── lambda.fa.{amb,ann,bwt,pac,sa,fai}
    ├── bam/
    │   ├── DTP1-2.sorted.bam
    │   ├── DTP1-2.sorted.bam.bai
    │   ├── DTP2-2.sorted.bam
    │   └── DTP2-2.sorted.bam.bai
    ├── stats/
    │   ├── DTP1-2.flagstat.txt
    │   └── DTP2-2.flagstat.txt
    └── coverage/
        ├── DTP1-2.samtools-coverage.tsv
        ├── DTP1-2.bamtocov.bed
        ├── DTP2-2.samtools-coverage.tsv
        └── DTP2-2.bamtocov.bed

(plus the files that were created in the genome/ directory by datasets: README.md, md5sum.txt, lambda.zip and ncbi_dataset/)

BAM files are big: if you need to free some space, you can delete genome/lambda.zip and the ncbi_dataset/ directory (once you have lambda.fa).

8. Challenge: minimap2

We installed minimap2 too. It is mainly used for long reads (Nanopore), but it can map short reads with the preset -ax sr. Map the same reads of DTP1-2 with it, writing a different BAM file, and compare the number of mapped reads with bwa:

SAMPLE=DTP1-2
minimap2 -t 4 -ax sr \
  $WORK/mapping/genome/lambda.fa \
  $(ls $READS/lab/PID-2959-${SAMPLE}_S*_R1_001.fastq.gz) \
  $(ls $READS/lab/PID-2959-${SAMPLE}_S*_R2_001.fastq.gz) \
  | samtools sort -@ 4 -o $WORK/mapping/bam/${SAMPLE}.minimap2.sorted.bam -

samtools flagstat $WORK/mapping/bam/${SAMPLE}.minimap2.sorted.bam

Which tool maps more reads? Why do you think the results are different?