NGS preprocessing and quality control

Overview

Objectives

• Inspect basic read statistics and quality profiles with FastQC
• Identify an overrepresented sequence and remove it with Cutadapt
• Trim low-quality bases with Trimmomatic and evaluate the effect with FastQC
• Handle paired-end reads as a Galaxy paired collection
• Compare Illumina, 454 and Ion Torrent read-length and quality profiles

Time: about 3 hours

Introduction

Today we will preprocess several types of NGS (Next Generation Sequencing) data. Preprocessing is the first step in almost any downstream analysis and ensures that we are working with high-quality sequencing reads.

Raw sequencing data often contain issues such as low-quality bases, sequencing adapters, or other technical artifacts. Identifying and correcting these problems improves the reliability of later analyses such as genome assembly, read alignment, and variant detection.

In this exercise we will work with datasets produced using different sequencing technologies and sequencing strategies.

Datasets used in this exercise:

1. Pseudomonas aeruginosa — single-end Illumina reads
2. Human — Illumina paired-end reads
3. Escherichia coli — 454 and Ion Torrent reads

P. aeruginosa single-end Illumina reads

Let us start with sequencing reads from a study of Pseudomonas aeruginosa, an opportunistic pathogen that can live in the environment and may infect humans. In humans it often colonizes the lungs, where it can form biofilms.

The reads come from a strain that infected approximately 40 patients in Denmark over the last 30 years. During this time the strain adapted to the human lung environment.

In a later alignment exercise (not part of today's exercise) we will compare these reads to the PAO1 reference genome. The PAO1 strain was isolated in 1955 in Australia from a wound and likely originated from the environment.

By comparing the genomes we can investigate which genes may have been lost during adaptation to the human host, assuming that the PAO1 genome represents a relatively close ancestor.

Before performing any downstream analysis we must first examine the quality of the sequencing data.

Step 1 — Get the data into Galaxy

All the data for today's exercise are collected in one Galaxy history. You will import that history once and use it for all three parts of the exercise.

Hands-on: Import the data

1. Log in to Galaxy Europe (create a free account first if you do not have one). Log in before you open the data link: otherwise your copy is not saved in your account.
2. Open the shared data history: shared data history
3. Click Import this history (top right, above the history).
4. In the Copying History dialog:
   • "Enter a title for the new history": replace the suggested title with Preprocessing
   • keep Copy only the active, non-deleted datasets. selected
   • click Copy History.
5. A blue message confirms "History imported and is now your active history", but the page still shows the shared history. Click the Galaxy logo at the top left (Home) to return to the main page. Your copy Preprocessing is now shown in the history panel on the right. If the panel shows another history, open Histories in the left-hand activity bar and switch to Preprocessing.
6. Check that the history contains the five files Paeruginosa.fastq.gz, HG00418_A_1.fastq.gz, HG00418_A_2.fastq.gz, GOS1.fastq.gz and iontorrent.fastq.gz, and the paired collection HG00418 ("a pair with 2 fastqsanger.gz datasets"). The history also shows 2 hidden datasets: these are the collection's two files and can be ignored.
7. Click on Paeruginosa.fastq.gz to expand it. Check that its datatype is fastqsanger.gz. A preview of the first reads is shown below the dataset details.

Tip: Only if the shared history does not work: upload the files yourself

If you cannot use the shared history, download the five files listed on the Data page to your computer and upload them to a new Galaxy history:

1. Click Upload in the left-hand activity bar.
2. Click Choose local file and select the files.
3. Set Type (set all) to fastqsanger.gz.
4. Click Start and wait until all datasets turn green.
5. Build the HG00418 paired collection as described in the paired-end part of this exercise.

Note: What does the datatype fastqsanger.gz mean?

Galaxy uses the datatype to decide which tools can read a dataset. fastqsanger means a FASTQ file whose quality scores use the "Sanger" encoding, and .gz means it is gzip-compressed. Trimmomatic only accepts inputs whose datatype states the quality encoding, so the datatype must be correct before trimming.

Step 2 — Quality control with FastQC

Quality control is an iterative process

Quality control of sequencing data is rarely a single step. Instead, it is typically an iterative process:

• Inspect the raw sequencing reads using FastQC
• Identify potential problems (low-quality bases, adapters, sequence biases)
• Apply trimming or filtering to correct these problems
• Run FastQC again to evaluate whether the quality has improved

This process may be repeated several times until the read quality is acceptable for downstream analyses such as genome assembly, alignment, or variant detection.

Hands-on: Run FastQC on the raw reads

1. Run FastQC (open tool) with the following parameters:
   • "Raw read data from your current history": Paeruginosa.fastq.gz (Galaxy pre-selects the newest dataset, so change it)
   • Leave all other parameters at their defaults.
   • Galaxy Europe shows a banner recommending the tool Falco instead. Ignore it: this exercise uses FastQC.
2. Click Run Tool.
3. FastQC produces two datasets: a Webpage (the report) and RawData (the same results as text).
4. When the job has finished (green), click the eye icon on the Webpage dataset to view the report.

Tip: If a job turns red

A red dataset means the job failed. See If a job fails on the course's Galaxy page. Temporary server problems are fixed by clicking Run Job Again and then Run Tool.

Tip: Reading about the tool

While the job is running, scroll down in the FastQC tool form to read the tool's help text.

Basic statistics

The first table in the FastQC report, Basic Statistics, summarises the whole file.

Questions

Q1. How many reads are there?

Q2. How long are the reads?

Q3. What is the GC percentage?

Q4. Which quality score encoding has been used? What is the corresponding score offset? (This information is needed when trimming reads or performing alignments.)

Quality plots and overrepresented sequences

The FastQC report contains several plots summarizing the quality of the sequencing reads. Look at each figure and check whether any warnings or failures are reported (the icons in the Summary list at the top of the report: a green tick means pass, an orange exclamation mark means warning, and a red cross means failure).

The most important plot is Per base sequence quality. This plot shows the quality score for each base position as you move from the 5′ end to the 3′ end of the reads. The yellow boxes show the range where the middle 50 % of quality values lie, and the background colours mark good (green), reasonable (orange) and poor (red, below 20) quality.

For Illumina data it is common to see a drop in quality toward the 3′ end of the reads. These low-quality bases should be removed because they may negatively affect downstream analyses such as genome assembly or SNP detection.

Questions

Q5. Look at the Per base sequence quality plot.

• How does the quality change from the 5′ end to the 3′ end of the reads?
• At approximately which position do the yellow boxes (the middle 50 % of quality values) start to reach into the red zone (below quality 20)?
• What happens at the last few positions of the reads?

Also examine the Overrepresented sequences table. These sequences are often adapters or PCR primers that remain in the reads after sequencing.

Adapters must be removed because they can:

• cause incorrect overlaps during de novo assembly
• interfere with read alignment
• introduce artificial sequence similarity between reads

In this dataset you should find an overrepresented sequence that FastQC identifies as an Illumina PCR primer (see the Possible Source column).

Record the sequence exactly as FastQC reports it. You will need it in the next step.

Questions

Q6. Which overrepresented sequence does FastQC report? What does the Possible Source column say, and how many times does the sequence occur?

Tip: Copying the sequence

Select the sequence in the FastQC report and copy it (Ctrl+C / Cmd+C). Check that you have copied only the sequence, with no spaces. Cutadapt accepts upper-case sequences only.

Questions: Predict before you trim

Q6b. Based on your FastQC report, which problem(s) should be addressed before any downstream analysis? For each problem, name the FastQC panel that shows it. What do you think the overrepresented sequence is, and in which part of the reads would you expect to find it?

Step 3 — Remove adapters using Cutadapt

Next we will remove the primer sequence from the reads using the program Cutadapt. We will search for it in two ways:

• as a 3′ adapter, which can occur anywhere in the read, including only partially at the 3′ end. Everything from the adapter to the end of the read is removed.
• as an anchored 5′ adapter, which must occur in full at the very beginning of the read. In Cutadapt a 5′ adapter is anchored by writing a ^ in front of the sequence.

Hands-on: Run Cutadapt

1. Run Cutadapt (open tool) with the following parameters:
   • "Single-end or Paired-end reads?": Single-end
     • "FASTQ/A file": Paeruginosa.fastq.gz
     • In "Read 1 Adapters":
       • "Insert 3' (End) Adapters"
         • "Source": Enter custom sequence
           • "Custom 3' adapter sequence": the sequence you recorded from FastQC. The box already contains an example value (AATTGGCC): delete it first, so that the box contains only your sequence.
       • "Insert 5' (Front) Adapters"
         • "Source": Enter custom sequence
           • "Custom 5' adapter sequence": ^ followed immediately by the same sequence (this makes the 5′ adapter anchored). Again, delete the example value first.
   • In "Adapter Handling Options": keep the defaults ("Maximum error rate": 0.1, "Minimum overlap length": 3).
   • In "Read Filtering Options" (click the section title to expand it):
     • "Minimum length (R1)": 50
   • "Additional outputs to generate": tick Report
2. Click Run Tool.
3. Rename the output Read 1 Output to Paeruginosa.cut: click the dataset, click the pencil icon, change Name, and click Save.

After trimming, reads shorter than 50 nt are discarded ("Minimum length (R1)"). This and the other Cutadapt settings can be changed in the tool form. More information about adapter types can be found in the Cutadapt documentation — Adapter types.

Inspect the Cutadapt report

When the job has finished, click the eye icon on the Report dataset. The statistics are shown at the top of the report (=== Summary ===), followed by one section per adapter.

Questions

Q7. How many reads contained the adapter (Reads with adapters)? Compare this with the number of occurrences FastQC reported. If the numbers differ, suggest why. (Hint: look at the length column in the table "Overview of removed sequences" — how long are most of the removed sequences?)

Q8. How many reads were removed because they were too short?

Q7b. Did the Cutadapt result support your prediction in Q6b? What do these numbers tell you about how much of the original dataset was affected by adapter/primer sequence?

Step 4 — Trim low-quality bases with Trimmomatic

Next we will remove low-quality bases from the ends of the reads using Trimmomatic. Trimmomatic applies its trimming operations in the order you list them, so add them in the order given below.

Questions: Predict before you trim

Q8b. Look again at the Per base sequence quality plot of the raw reads (Q5). Which part of the reads do you expect quality trimming to affect most? What do you expect to happen to the read lengths, and could some reads disappear completely?

Hands-on: Run Trimmomatic

1. Run Trimmomatic (open tool) with the following parameters:
   • "Single-end or paired-end reads?": Single-end
     • "Input FASTQ file": Paeruginosa.cut
   • "Perform initial ILLUMINACLIP step?": no
   • In "Trimmomatic Operation", set the first operation and use Insert Trimmomatic Operation to add the others:
     1. "Select Trimmomatic operation to perform": Cut bases off the start of a read, if below a threshold quality (LEADING)
        • "Minimum quality required to keep a base": 20
     2. Cut bases off the end of a read, if below a threshold quality (TRAILING)
        • "Minimum quality required to keep a base": 20
     3. Sliding window trimming (SLIDINGWINDOW)
        • "Number of bases to average across": 4
        • "Average quality required": 20
     4. Drop reads below a specified length (MINLEN)
        • "Minimum length of reads to be kept": 36
   • "Output trimmomatic log messages?": Yes
2. Click Run Tool.
3. Rename the trimmed reads (the output named Trimmomatic on Paeruginosa.cut, not the log file) to Paeruginosa.cut.trim.

The trimming statistics are in the log file output. Click its eye icon and find the line starting with Input Reads:.

Questions

Q9. How many reads survived the trimming? How many were dropped?

Q9b. Which of the four operations is the one designed to remove whole reads? Using your answer to Q8b, explain why so many more reads were lost here than in the Cutadapt step.

Step 5 — Compare FastQC reports

After trimming, examine how the dataset has changed.

Hands-on: Run FastQC on the trimmed data

1. Run FastQC (open tool):
   • "Raw read data from your current history": select both Paeruginosa.cut and Paeruginosa.cut.trim. Click the Multiple datasets button, then click the datasets in the list. Galaxy may already have selected one dataset for you, and clicking a selected dataset removes it again: check that the box shows exactly these two datasets (remove any other with its ×). Galaxy runs FastQC once for each dataset.
2. Open the three FastQC Webpage reports (raw, Paeruginosa.cut and Paeruginosa.cut.trim). Galaxy names each report after the number of its input in your history: for example, FastQC on dataset 1: Webpage is the report for dataset 1, Paeruginosa.fastq.gz.

Tip: Viewing reports side by side

Click the grid icon in the top bar (Enable/Disable Window Manager). Each report you open with the eye icon then appears in its own window, which you can move and resize to arrange the three reports next to each other. When you have finished, close the windows (×) and click the grid icon again to switch the Window Manager off; otherwise the windows cover the tool forms. Alternatively, open each report in its own browser tab: Ctrl-click (Cmd-click on a Mac) the eye icon.

Questions

Q10. Compare the three FastQC reports side by side: before (raw reads) and after each preprocessing step (adapter-trimmed and quality-trimmed reads). Look at the plots themselves, not only the summary icons.

1. How did the Per base sequence quality plot change after adapter trimming, and after quality trimming?
2. Is the overrepresented sequence still reported after adapter trimming?
3. Did the Per base sequence content plot change after trimming?
4. Which part of the reads still shows an unusual pattern in Per base sequence content? Suggest an additional trimming step, and which tool option could perform it.
5. Before vs after: did preprocessing solve the problems you identified in Q6b? For each step, say what improved, what did not materially change, and whether the step had its intended effect.

Human Illumina paired-end reads

Next we will examine paired-end Illumina reads. These reads originate from a ~40× whole-genome sequencing (WGS) dataset of an Asian individual.

To keep the exercise manageable, the dataset has been filtered so that it only contains reads mapping to chromosome 21.

Each DNA fragment was sequenced from both ends. The file HG00418_A_1.fastq.gz contains the forward (left) read of each pair, and HG00418_A_2.fastq.gz contains the reverse (right) read. The n-th read in one file is the mate of the n-th read in the other file. In Galaxy we keep the two files together as a paired collection, so that tools always receive the forward and reverse reads together and in the right order.

Step 1 — The paired collection

Hands-on: Check the paired collection

1. In your history, click the HG00418 collection. It contains two datasets, shown as forward (the reads from HG00418_A_1.fastq.gz) and reverse (the reads from HG00418_A_2.fastq.gz).
2. Click « History: Preprocessing at the top of the history panel to return to the full history.

Tip: If you uploaded the files yourself: build the paired collection

1. At the top of the history, click the Select Items (checkbox) icon.
2. Tick HG00418_A_1.fastq.gz and HG00418_A_2.fastq.gz.
3. Open the selection dropdown ("2 of … selected") and choose Advanced Build List.
4. Under What are you building? choose List of Paired Datasets and click Next.
5. Under Auto Pairing, check that the two files are paired with HG00418_A_1.fastq.gz as the forward read and HG00418_A_2.fastq.gz as the reverse read, and click Next.
6. Name the collection HG00418 and click Build.

This creates a list containing one pair. You can use it everywhere this exercise says HG00418; Galaxy then runs each tool on the pair inside the list.

Step 2 — Run FastQC

Hands-on: Run FastQC on the paired collection

1. Run FastQC (open tool):
   • "Raw read data from your current history": click the Dataset collection button and select HG00418.
2. Galaxy runs FastQC separately on the forward and on the reverse reads. The Webpage output is a collection with one report for each.

Questions

Q11. Use Basic Statistics in the two reports.

1. How many reads are in the forward file, and how many in the reverse file?
2. How many read pairs is that, and how many individual reads (mates) in total?
3. How long are the reads?

Q12. Compare the Per base sequence quality plots of the forward and reverse reads. How does the quality change toward the 3′ end, and do the two mates differ?

Q13. Does FastQC report any adapter contamination (Adapter Content) or overrepresented sequences in this dataset? What does that mean for how we should trim these reads?

Step 3 — Trim reads using Trimmomatic

Based on your answers above, we will now trim the paired-end reads with Trimmomatic, using the same quality settings as in the previous exercise.

Hands-on: Run Trimmomatic on the paired collection

1. Run Trimmomatic (open tool) with the following parameters:
   • "Single-end or paired-end reads?": Paired-end (as collection)
     • "Select FASTQ dataset collection with R1/R2 pair": HG00418
   • "Perform initial ILLUMINACLIP step?": no
   • "Trimmomatic Operation": the same four operations, in the same order, as for P. aeruginosa:
     1. LEADING: 20
     2. TRAILING: 20
     3. SLIDINGWINDOW: "Number of bases to average across" 4, "Average quality required" 20
     4. MINLEN: 36
   • "Output trimmomatic log messages?": Yes
2. Click Run Tool.
3. Trimmomatic produces two collections: paired (forward + reverse reads whose mate also survived) and unpaired (reads whose mate was removed), and a log file.

Open the log file and find the line starting with Input Read Pairs:.

Questions

Q14. Record the trimming statistics:

• Input read pairs
• Both surviving
• Forward only surviving
• Reverse only surviving
• Dropped

Q14b. Why can paired-end trimming produce these four categories? In this dataset, more forward only than reverse only reads survived. Which of your FastQC observations (Q12) explains that?

Interpretation

Unlike the previous dataset, this dataset contains paired-end reads. Each DNA fragment was sequenced from both ends, producing a forward and a reverse read.

During trimming, some reads may fail the quality filters while their partner read remains usable.

This is why the output contains categories such as:

• Both surviving — both reads in the pair passed filtering (Galaxy output: the paired collection)
• Forward only surviving — only the forward read passed (Galaxy output: forward in the unpaired collection)
• Reverse only surviving — only the reverse read passed (Galaxy output: reverse in the unpaired collection)
• Dropped — both reads failed filtering

Reads where only one partner survives are no longer paired and are therefore treated as single-end reads in downstream analyses.

Questions

Q15. How many reads in total ended up in the unpaired collection? Could they be used together with the paired reads in a paired-end alignment?

E. coli 454 and Ion Torrent data

In this final exercise we will examine sequencing data generated using two older sequencing technologies: 454 pyrosequencing and Ion Torrent.

The reads originate from the German E. coli outbreak in 2011 associated with contaminated cucumbers and sprouts.

Note that these datasets were not generated using Illumina sequencing. In Galaxy they still have the datatype fastqsanger.gz, because the datatype describes the file format and quality encoding, not the sequencing technology. Keep track of which technology produced which file by giving the datasets clear names.

Step 1 — The data

Hands-on: Name the datasets

1. Find GOS1.fastq.gz (the 454 dataset) and iontorrent.fastq.gz (the Ion Torrent dataset) in your history.
2. Optional: rename them to GOS1 (454) and iontorrent (Ion Torrent) so that the FastQC reports are easy to tell apart.

Step 2 — Run FastQC

Hands-on: Run FastQC on both datasets

1. Run FastQC (open tool):
   • "Raw read data from your current history": select both GOS1.fastq.gz and iontorrent.fastq.gz (Multiple datasets). As before, check that the box shows exactly these two datasets.
2. Open both Webpage reports.

Pay particular attention to Basic Statistics, the Per base sequence quality plot and the Sequence Length Distribution.

Questions

Q16. For each dataset, record from Basic Statistics:

• number of reads
• read length (range)
• GC percentage

Q17. Open the Per base sequence quality, Per sequence quality scores and Sequence Length Distribution plots of both reports side by side. What differences do you observe in read length, in quality, and in how quality deteriorates along the reads?

Q18. Compare all three technologies you have seen today (Illumina, 454, Ion Torrent). Which produce reads of a fixed length? Summarise, in one sentence each, the characteristic read length and quality profile of each technology. Based on the plots, which dataset would you expect to lose the most reads if trimmed with the same settings as the Illumina data (average quality 20), and why? (You can test your prediction in the optional exploration, Q19.)

Optional exploration

Try trimming the 454 and Ion Torrent reads to remove low-quality bases with Trimmomatic (open tool) ("Single-end", the same four operations as before), and then run FastQC again to evaluate the effect of trimming.

Questions

Q19. (Optional) How many reads survive trimming for each technology? Explain what you observe.

Putting it together

Quality control is a loop: observe → predict → process → observe again → compare → explain.

Questions

Q20. You receive a new sequencing dataset tomorrow. Before any downstream analysis:

1. What would you inspect first, and in which FastQC panels?
2. How would you decide whether trimming is necessary, and which kind?
3. What would you check after preprocessing?

Exercise completed

Congratulations — you have now completed the data preprocessing exercises.

You have explored sequencing data generated using multiple technologies and learned how to:

• inspect sequencing read statistics
• evaluate read quality using FastQC
• remove adapters and primers
• trim low-quality bases
• compare different sequencing technologies