Overview

Time estimation: HM

Version: main

Last update: 2026-05-29

Questions:

• What is quality processing and why it is important?

• What are the neccessary steps to do quality processing on 16S amplicon data?

Objectives:

• You will be able to distinguish between different types of errors and bias in amplicon data.

• You will apply different tools to address the different error types on the data.

• You will compare the effect of quality processing with unporcessed data. -

The Tutorial Data Set

The first thing you need to do is to connect to your virtual machine with the X2Go Client.
If you are working with your laptop and haven’t installed it yet, you can get it here:

https://wiki.x2go.org/doku.php/download:start

1. Enter the IP of your virtual machine, the port, the username ubuntu, and select your SSH key.
2. When you have successfully connected to your machine, open a terminal.

Give the attached volume to the `ubuntu` user

sudo chown ubuntu:ubuntu /mnt/volume/

Create a link to the mounted volume in your home directory

ln -s /mnt/volume/ workdir

Download the tutorial dataset (and pre‑computed results)

cd ~/workdir
wget https://openstack.cebitec.uni-bielefeld.de:8080/swift/v1/mgcourse_data/16S-data.tgz

Unpack the archive

tar -xzvf 16S-data.tgz
ln -s 16S-data/raw/ 16Sdata

Clean up

rm 16S-data.tgz

Inspect the data directory

(mgcourse) ubuntu@sebmain-939e0:~/workdir$ ls -la 16Sdata/
total 179704
drwxrwxr-x 2 ubuntu ubuntu     4096 Oct  8 12:12 .
drwxr-xr-x 4 ubuntu ubuntu     4096 Nov 20  2015 ..
-rw-rw-r-- 1 ubuntu ubuntu  5712157 Nov 20  2015 BGA1_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  5712157 Nov 20  2015 BGA1_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 21370972 Nov 20  2015 BGA1_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 21370972 Nov 20  2015 BGA1_2_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu  8910118 Nov 20  2015 BGA2_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  8910118 Nov 20  2015 BGA2_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 11004978 Nov 20  2015 BGA2_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 11004978 Nov 20  2015 BGA2_2_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 19342508 Nov 20  2015 BGA3_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 19342508 Nov 20  2015 BGA3_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 12128141 Nov 20  2015 BGA3_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 12128141 Nov 20  2015 BGA3_2_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu  9294674 Nov 20  2015 BGA4_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  9294674 Nov 20  2015 BGA4_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu  4222729 Nov 20  2015 BGA4_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  4222729 Nov 20  2015 BGA4_2_R2.fastq
-rw-r--r-- 1 ubuntu ubuntu     1240 Nov 23  2015 Pipe.tsv
-rw-rw-r-- 1 ubuntu ubuntu       55 Nov 20  2015 Primers.txt
(mgcourse) ubuntu@sebmain-939e0:~/workdir$

Inspect the data directory

ls -la 16Sdata/

Code-out

total 179704
drwxrwxr-x 2 ubuntu ubuntu     4096 Oct  8 12:12 .
drwxr-xr-x 4 ubuntu ubuntu     4096 Nov 20  2015 ..
-rw-rw-r-- 1 ubuntu ubuntu  5712157 Nov 20  2015 BGA1_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  5712157 Nov 20  2015 BGA1_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 21370972 Nov 20  2015 BGA1_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 21370972 Nov 20  2015 BGA1_2_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu  8910118 Nov 20  2015 BGA2_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  8910118 Nov 20  2015 BGA2_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 11004978 Nov 20  2015 BGA2_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 11004978 Nov 20  2015 BGA2_2_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 19342508 Nov 20  2015 BGA3_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 19342508 Nov 20  2015 BGA3_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu 12128141 Nov 20  2015 BGA3_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu 12128141 Nov 20  2015 BGA3_2_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu  9294674 Nov 20  2015 BGA4_1_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  9294674 Nov 20  2015 BGA4_1_R2.fastq
-rw-rw-r-- 1 ubuntu ubuntu  4222729 Nov 20  2015 BGA4_2_R1.fastq
-rw-rw-r-- 1 ubuntu ubuntu  4222729 Nov 20  2015 BGA4_2_R2.fastq
-rw-r--r-- 1 ubuntu ubuntu     1240 Nov 23  2015 Pipe.tsv
-rw-rw-r-- 1 ubuntu ubuntu       55 Nov 20  2015 Primers.txt
(mgcourse) ubuntu@sebmain-939e0:~/workdir$

(Optional) Disable system beep sounds

xset -b

Activate Conda Environment

In order to follow the hands‑on session you need to activate a special Conda environment where all tools and software have been installed and are available:

conda activate mgcourse

Note: You need to run the command above in every new terminal you open.

FASTQ Format

Let’s first have a look at what a FASTQ file looks like and how its format is defined:

A FASTQ file contains four lines per sequence:

1. Sequence identifier, starting with @
2. Raw sequence letters
3. A separator line, beginning with + (can optionally be followed by the identifier)
4. Quality string, encoding per-base quality scores

Sequence Quality Scores

Sequence quality scores were first introduced by the phred base caller used for Sanger sequencing.

Quality scores (Q-scores) are a log transformation of the error probability for each base:

Quality value	Error probability
20	1/100
30	1/1000
40	1/10,000
50	1/100,000

• Quality score = log-transformed error probability

The higher the quality score, the lower the probability that the base call is incorrect.

References

Images: CUNY OpenLab, Next-gen sequencing

FastQC

FastQC provides an easy way to perform quality control checks on raw high-throughput sequencing data. Its modular set of analyses allows you to quickly determine if your data have problematic features that could affect downstream analysis.

Main functions of FastQC:

• Import data from BAM, SAM, or FASTQ files (any variant)
• Provide a quick overview to identify potential problems
• Summarize data with graphs and tables
• Export results into permanent HTML reports
• Allow offline/automated report generation (including CLI)

You can run FastQC interactively or via the command-line interface (CLI).

Running FastQC

To see all FastQC options, use:

fastqc --help

Common usage:

Code-in

fastqc seqfile1 seqfile2 .. seqfileN
fastqc [-o output_dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

Frequently used options:

Option	Description
-o, --outdir	Output directory (must exist)
--casava	Treat files as raw Casava output (groups samples)
--nano	Special support for Nanopore fast5 directories
--nofilter	With –casava, do not remove poor quality reads
--nogroup	Disable base grouping for long reads (be careful!)
-f, --format	Force a specific file format (fastq, bam, sam, etc)
-t, --threads	Number of threads/files processed in parallel
-c, --contaminants	File with named contaminants to check for
-a, --adapters	File with named adapter sequences to check for
-l, --limits	Custom file with module warn/error thresholds
-k, --kmers	Kmer length for kmer content module (default 7)
-q, --quiet	Suppress progress messages
-d, --dir	Temporary file directory for images

For more detailed information, see the FastQC documentation.

Evaluating Quality with FastQC

Let’s evaluate a set of FASTQ files using FastQC from the command line.

1. Create a results directory and run FastQC on your files.

   cd ~/workdir
   mkdir -p ~/workdir/fastqc  
   fastqc -t 14 -o ~/workdir/fastqc/ ~/workdir/16Sdata/*.fastq
   

2. View the reports in your web browser (replace *.html with your actual report name if needed):

   cd ~/workdir/fastqc
   firefox *.html
   

3. Browse through the FastQC HTML report. Together, we will inspect which modules flag your data as poor, warning, or pass.

FastQC Example Reports

See the FastQC homepage for example reports—some with notably bad quality!

References

• FastQC

Quality Treatment

• Reads contain errors (0.1–15%) and contamination
• Quality matters!
  • NGS high throughput = lots of data → more data, more (systematic) errors
  • 16S data > WGS read-based > WGS assembly-based (in increasing stringency)
  • Better data = lower computational cost
  • Decreases false positives
  • Always a trade-off (risk of false negatives)
• Quality-based filtering vs. error correction
• Many tools available for this step

Typical Workflow

For this exercise, we will use:

• Only 16S data
• Merge: FLASh
• Clip primers: cutadapt
• Trim quality: sickle
• Filter by length: ea-utils

Merge Reads

• Assembly of forward and reverse read pairs (if DNA fragment is shorter than twice the read length)
• Ungapped alignment with min overlap region (common to Illumina)
• Quality scores at merged positions are recalculated as absolute difference

Merge first pair of reads

Code-in

mkdir -p ~/workdir/flash
cd ~/workdir/flash
flash -r 300 ~/workdir/16Sdata/BGA1_1_R1.fastq ~/workdir/16Sdata/BGA1_1_R2.fastq -o BGA1_1

Example output:

[FLASH] Read combination statistics:
[FLASH]     Total pairs:      8495
[FLASH]     Combined pairs:   8388
[FLASH]     Uncombined pairs: 107
[FLASH]     Percent combined: 98.74%

The merged paired-end reads have been written to BGA1_1.extendedFrags.fastq

Merge all pairs in batch

Code-in

parallel "flash -r 300 ~/workdir/16Sdata/{}_R1.fastq ~/workdir/16Sdata/{}_R2.fastq -o {}" ::: {BGA1_2,BGA2_1,BGA2_2,BGA3_1,BGA3_2,BGA4_1,BGA4_2}

If you have information on the amplified fragment, you can adjust min/max overlap, fragment length, and SD as necessary.

Primer Clipping

First, identify the primer sequences used to amplify your region of interest:

Domain	Region	Sequence
Bacteria	V3F	CTACGGGNGGCWGCAG
Bacteria	V4R	GACTACHVGGGTATCTAATCC

Clip forward primer with cutadapt

Code-in

mkdir ~/workdir/cutadapt
cd ~/workdir/cutadapt
cutadapt -g ^CTACGGGNGGCWGCAG ../flash/BGA1_1.extendedFrags.fastq -e 0.2 -O 10 -o BGA1_1.trimmedf.fastq

• ^ anchors the primer at the 5’ end of the read
• -e 0.2 sets a max error rate of 20%
• -O 10 is a minimum overlap of 10 bases
• Cutadapt accepts wobble bases and can trim both primers and adapters
• Adjust stringency parameters as needed; always inspect output for too many/suspiciously trimmed reads!

Example cutadapt output:

=== Summary ===
Total reads processed:                   8,388
Reads with adapters:                     8,383 (99.9%)
Reads written (passing filters):         8,388 (100.0%)
Total basepairs processed:     3,808,321 bp
Total written (filtered):      3,665,848 bp (96.3%)
=== Adapter 1 ===
Sequence: CTACGGGNGGCWGCAG; Type: anchored 5'; Length: 16; Trimmed: 8383 times.
No. of allowed errors:
0-4 bp: 0; 5-9 bp: 1; 10-14 bp: 2; 15 bp: 3

Clip reverse primer (reverse complement!)

After merging, the reverse primer must be reverse complemented.

1. Reverse-complement the primer:

   cd ~/workdir
   echo -e ">primer\nGACTACHVGGGTATCTAATCC" > revprimer.fas
   revseq -sequence revprimer.fas -outseq revprimer_rc.fas
   cat revprimer_rc.fas

2. Use the reverse-complemented primer sequence with `cutadapt` at the 3' end:

   cd ~/workdir/cutadapt
   cutadapt -a GGATTAGATACCCBDGTAGTC$ BGA1_1.trimmedf.fastq -e 0.2 -O 10 -o BGA1_1.trimmedfr.fastq

Batch primer trimming

Code-in

cd ~/workdir/cutadapt
parallel "cutadapt -g ^CTACGGGNGGCWGCAG ../flash/{}.extendedFrags.fastq -e 0.2 -O 10 -o {}.trimmedf.fastq" ::: {BGA1_2,BGA2_1,BGA2_2,BGA3_1,BGA3_2,BGA4_1,BGA4_2}
parallel "cutadapt -a GGATTAGATACCCBDGTAGTC$ {}.trimmedf.fastq -e 0.2 -O 10 -o {}.trimmedfr.fastq" ::: {BGA1_2,BGA2_1,BGA2_2,BGA3_1,BGA3_2,BGA4_1,BGA4_2}

Quality Trimming

Reads with very low quality often contain many miscalled bases, which can increase false positives in downstream analyses. Trim low-quality bases from both ends using a sliding window approach with sickle.

Trim low-quality ends

Code-in

mkdir -p ~/workdir/sickle
cd ~/workdir/sickle
sickle se -f ../cutadapt/BGA1_1.trimmedfr.fastq -t sanger -q20 -o BGA1_1.clipped.fastq

• -q 20 sets min average quality score to 20
• -t sanger applies Phred+33 scale
• -n truncates at ambiguous bases (N)

Batch quality trimming

Code-in

cd ~/workdir/sickle
parallel "sickle se -f ../cutadapt/{}.trimmedfr.fastq -t sanger -q20 -o {}.clipped.fastq" ::: {BGA1_2,BGA2_1,BGA2_2,BGA3_1,BGA3_2,BGA4_1, BGA4_2}

Length Filtering

Filter out reads that are too short (or too long) for your expected amplicon size. First, use a Perl script to inspect read lengths:

Plot and inspect read lengths

Code-in

cd $CONDA_PREFIX/bin  
wget https://raw.githubusercontent.com/jueneman/16S-workshop-denbi/master/docs/qc/FastaStats.pl
chmod u+x FastaStats.pl

mkdir ~/workdir/length
cd  ~/workdir/length  
FastaStats.pl -q ../sickle/BGA1_1.clipped.fastq > BGA1_1.fastq.hist
head -n 10 BGA1_1.fastq.hist

Now, trim reads to a defined length window (here, using 1.5 × IQR):

Filter by length with ea-utils (fastq-mcf)

Code-in

fastq-mcf -0 -l 369 -L 461 n/a ../sickle/BGA1_1.clipped.fastq -o BGA1_1.fastq

• -l 369 = lower length threshold
• -L 461 = upper length threshold
• n/a = do not provide primer file (no adapter trimming here)

Batch length filtering

Code-in

parallel "FastaStats.pl -q ../sickle/{}.clipped.fastq > {}.fastq.hist" ::: {BGA1_2,BGA2_1,BGA2_2,BGA3_1,BGA3_2,BGA4_1,BGA4_2}

grep IQR *.hist
fastq-mcf -0 -l 369 -L 461 n/a ../sickle/BGA1_2.clipped.fastq -o BGA1_2.fastq
parallel "fastq-mcf -0 -l 372 -L 460 n/a ../sickle/{}.clipped.fastq -o {}.fastq" ::: {BGA2_1,BGA2_2}
parallel "fastq-mcf -0 -l 364 -L 464 n/a ../sickle/{}.clipped.fastq -o {}.fastq" ::: {BGA3_1,BGA3_2}
parallel "fastq-mcf -0 -l 377 -L 449 n/a ../sickle/{}.clipped.fastq -o {}.fastq" ::: {BGA4_1,BGA4_2}

FastQC - Revisited

Now inspect how your quality treatment has improved the read quality:

Run FastQC on filtered files

Code-in

cd ~/workdir/fastqc
fastqc -t 14 -o ~/workdir/fastqc/ ~/workdir/length/*.fastq

As one final step, group all high-quality files together:

Organize your high-quality reads

Code-in

cd ~/workdir
mkdir HQ
cp length/*.fastq HQ/

Final Remarks

• Know your data: library prep, expected fragment/read length, possible adapters/primers
• Consider the sequencer (Illumina vs. Ion Torrent, etc.)
• Carefully inspect results at each step
• Try different strategies (conservative vs. loose parameters)
• Adapt the QC workflow to your research question (16S vs. read-based shotgun, etc.)
• This is only one workflow; consider other tools or workflow order. Always examine your raw data!

References

• FastQC: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
• Sickle: https://github.com/najoshi/sickle
• cutadapt: https://github.com/marcelm/cutadapt
• FLASh: http://ccb.jhu.edu/software/FLASH/
• ea-utils: https://github.com/ExpressionAnalysis/ea-utils

Key Points