Today we practise slowly on one pair of samples. For every step: inspect the input, run the command, then inspect the output.

• Raw reads: /shared/team/2026_training/week2/reads
• Your working directory: /shared/team/2026_training/week2/{firstname}
• Work only inside your own directory. Never write in reads/.

1. Your directory and notebook

Replace {firstname} with your directory name:

cd /shared/team/2026_training/week2/{firstname}
pwd
ls -lh
nano notebook.md

Write a header in the notebook (Ctrl+O, Enter, Ctrl+X to save):

# Day 4 notebook

Name: {firstname}

Your notebook is also your command file. Before running a command, write it there, check the paths, then copy it into the terminal. Do not use variables: write the full names, so the command is easy to read and correct. After each step, add what you observed.

cat notebook.md

2. Environment

The week2 environment contains the tools. Create it if you do not have it yet:

mamba create -n week2 -c conda-forge -c bioconda seqfu fastp kraken2 multiqc kraut
conda activate week2
seqfu version
fastp --version
kraken2 --version

Your prompt must start with (week2). You need conda activate week2 in every new terminal.

3. Find your four files

Find your sample number N:

Your sample was sequenced in pools DTP1 and DTP2. Below we use 4 as an example: replace every 4 with your number. The file names also contain a sequencing number (_S9_, _S17_…) that we do not know, so we use *.

ls /shared/team/2026_training/week2/reads/lab/*DTP1-4*
ls /shared/team/2026_training/week2/reads/lab/*DTP2-4*

You should see four files (R1 and R2 for each pool). Write their names in the notebook. Now inspect one raw file:

seqfu less /shared/team/2026_training/week2/reads/lab/PID-2959-DTP1-4_S*_R1_001.fastq.gz

Space moves, / searches, Q quits. Look at a few reads: name, sequence and quality. Do the same for R2.

Notebook: the four file names, and what a FASTQ record looks like.

4. Read statistics with seqfu stats

mkdir -p seqfu-stats

seqfu stats --basename /shared/team/2026_training/week2/reads/lab/PID-2959-DTP1-4_S*_R*_001.fastq.gz /shared/team/2026_training/week2/reads/lab/PID-2959-DTP2-4_S*_R*_001.fastq.gz > seqfu-stats/illumina.tsv

--basename prints only the file name. > writes the output to the file (and replaces it if it already exists). Inspect the result:

ls -lh seqfu-stats
cat seqfu-stats/illumina.tsv

The table has a header and one line per file. Columns: File, #Seq, Total bp, Avg, N50, N75, N90, auN, Min, Max.

Notebook:

1. How many reads are in each file? Do R1 and R2 have the same number?
2. Which pool has more reads? What is the typical read length?

5. Filtering with fastp

fastp removes adapters and low-quality reads. With no extra options it uses its default filters, which is enough today. Create the output directories:

mkdir -p filtered-reads/html-reports filtered-reads/json

Run it on DTP1:

fastp \
  --in1 /shared/team/2026_training/week2/reads/lab/PID-2959-DTP1-4_S*_R1_001.fastq.gz \
  --in2 /shared/team/2026_training/week2/reads/lab/PID-2959-DTP1-4_S*_R2_001.fastq.gz \
  --out1 filtered-reads/DTP1-4_R1.fastq.gz \
  --out2 filtered-reads/DTP1-4_R2.fastq.gz \
  --html filtered-reads/html-reports/DTP1-4.fastp.html \
  --json filtered-reads/json/DTP1-4.fastp.json

Now copy the command from your notebook and change every DTP1 to DTP2. Inspect the output:

ls -lh filtered-reads filtered-reads/html-reports filtered-reads/json
seqfu stats --basename filtered-reads/DTP1-4_R1.fastq.gz filtered-reads/DTP1-4_R2.fastq.gz
seqfu less filtered-reads/DTP1-4_R1.fastq.gz

Open filtered-reads/html-reports/DTP1-4.fastp.html in JupyterLab and compare before and after filtering.

Notebook:

1. How many reads passed? What percentage was discarded?
2. Were adapters removed?
3. Do R1 and R2 still have the same number of reads?

6. kraken2: what is in the sample?

Kraken2 compares the reads with a shared database of reference genomes.

mkdir -p kraken-reports

kraken2 \
  --db /shared/public/db/kraken2/pluspfp_16gb/2026-06-26/ \
  --paired \
  --threads 4 \
  --memory-mapping \
  --output /dev/null \
  --report kraken-reports/DTP1-4.tsv \
  filtered-reads/DTP1-4_R1.fastq.gz \
  filtered-reads/DTP1-4_R2.fastq.gz

Repeat for DTP2, changing DTP1 to DTP2. This can take a few minutes. --paired says that R1 and R2 belong together, --report is the summary table, --output /dev/null discards the per-read classification, and --memory-mapping shares the database between users.

Inspect the report:

ls -lh kraken-reports
head -n 20 kraken-reports/DTP1-4.tsv

It is sorted by abundance. Column 1 is the percentage of reads, column 2 the number of reads, column 4 the rank, and the last column the name. Rank U means unclassified, S means species.

Notebook:

1. What percentage of reads is unclassified in each pool?
2. What is the most abundant species? Do the two pools agree?
3. Which organism do you think you sequenced?

7. Final check

{firstname}/
├── notebook.md
├── seqfu-stats/illumina.tsv
├── filtered-reads/        (4 FASTQ, 2 HTML, 2 JSON)
└── kraken-reports/        (DTP1-N.tsv and DTP2-N.tsv)

find . -type f | sort

Checklist:

• notebook.md has the commands and your answers
• seqfu-stats/illumina.tsv has 5 lines (header + 4 files)
• filtered-reads/ has 4 FASTQ, 2 HTML and 2 JSON files
• kraken-reports/ has DTP1-N.tsv and DTP2-N.tsv

When you finish, let us know and keep your notebook open: we will use these files in the next steps.