Assessing Read Quality, Trimming and Filtering
Overview
Teaching: 45 min
Exercises: 45 minQuestions
How can I describe the quality of my data?
How can we get rid of sequence data that doesn’t meet our quality standards?
How do these methods differ when looking at Nanopore data?
Objectives
Interpret a FastQC plot summarizing per-base quality across all reads.
Interpret the NanoPlot output summarizing a Nanopore sequencing run
Filter Nanopore reads based on quality using the command line tool SeqKit
Quality control
Before assembling our metagenome from the the short-read Illumina sequences and the long-read Nanopore sequences, we need to apply quality control to both. The two types of sequence data require different QC methods. We will use:
- FastQC to examine the quality of the short-read Illumina data
- NanoPlot to examine the quality of the long-read Nanopore data and Seqkit to trim and filter them.
Illumina Quality control using FastQC
Reminder of the FASTQ format
In the FASTQ file format, each ‘read’ (i.e. sequence) is composed of four lines:
Line Description 1 Always begins with ‘@’ and gives the sequence identifier and an optional description 2 The actual DNA sequence 3 Always begins with a ‘+’ and sometimes the same info in line 1 4 Has a string of characters which represent the PHRED quality score for each of the bases in line 2; must have same number of characters as line 2
We can examine the first read in the FASTQ file using head to print the first four lines.
Move into the directory containing the illumina FASTQ files using cd:
cd cs_course/data/illumina_fastq/
Now use head to view the first 4 lines of the ERR2935805.fastq file
head -n 4 ERR2935805.fastq
@ERR2935805.1 HWI-C00124:284:HWTT2BCXY:1:1101:1247:2214 length=202
GATGGCGATAGAAGTCAAGTCTTTATTTTATGAAACCGCCATCATTAGTAGTATTTTATTTGGGCTCCCTTTTATAGGGACGGATATTTATGAGAATCAGCNAAAAAATCTACNCCTTCCTGAAANNANNAACNNNCAGGGTCTGACGATTTTCCTGCTGGGGTGGGAAATTGCCAGATAAAACAATATTGTGATTATCTCT
+ERR2935805.1 HWI-C00124:284:HWTT2BCXY:1:1101:1247:2214 length=202
GAGGGGIIIIIIGIIGIGIIGIGIGIIIIGIIGGGIGGGGGGGIIGIIIIIIIGGIGGIIIIGGGGGGIIGIIIGGGIIGGGGAGGGAGGGIAGGGGGGGA#<GGAIIIIIGI#<<GGGIGGGGA##<##<<G###<<GAGGIGGIIIIIIIIIGGIIIIIIIIIIGIIGGGGGAGGGGIIIGGIIGIGIGIGGIGGIGGG.
The quality score of this read is on line 4.
GAGGGGIIIIIIGIIGIGIIGIGIGIIIIGIIGGGIGGGGGGGIIGIIIIIIIGGIGGIIIIGGGGGGIIGIIIGGGIIGGGGAGGGAGGGIAGGGGGGGA#<GGAIIIIIGI#<<GGGIGGGGA##<##<<G###<<GAGGIGGIIIIIIIIIGGIIIIIIIIIIGIIGGGGGAGGGGIIIGGIIGIGIGIGGIGGIGGG.
PHRED score reminder
Quality encoding: !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJ | | | | | Quality score: 01........11........21........31........41Quality is interpreted as the probability of an incorrect base call. To make it possible to line up each individual nucleotide with its quality score, the numerical score is encoded by a single character. The quality score represents the probability that the corresponding nucleotide call is incorrect. It is a logarithmic scale so a quality score of 10 reflects a base call accuracy of 90%, but a quality score of 20 reflects a base call accuracy of 99%.
The PHRED quality scores for the majority of the bases in this read are between 31-41.
Rather than assessing every read in the raw data by hand we can use FastQC to visualise the quality of the whole sequencing file.
First, we are going to organise our analysis by creating a directory to contain the output of all of the analysis we generate in this course.
The mkdir command can be used to make a new directory. Using the -p flag for mkdir allows it to create a new directory, even if one of the parent directories doesn’t already exist. It also suppresses errors if the directory already exists, without overwriting that directory.
Return to your home directory (/home/csuser)
cd
Create the directories analysis/qc/illumina_qc inside cs_course
mkdir -p cs_course/analysis/qc/illumina_qc
You might want to use ls to check those nested directories have been made.
Now we have created the directories we are ready to start the quality control of the Illumina data.
FastQC has been installed on your instance and we can run it with the -h flag to display the help documentation and remind ourselves how to use it and all the parameters available:
$ fastqc -h
FastQC seq help documentation
We need to use just an input file and the -o flag to give the output directory: fastqc _inputfile_ -o _outputdirectory_
Navigate to your qc/ directory,
cd ~/cs_course/analysis/qc/
As we are using only one FASTQ file we can specify fastqc and then the location of the FASTQ file we want to analyse, and the illumina_qc output directory:
fastqc ~/cs_course/data/illumina_fastq/ERR2935805.fastq -o illumina_qc/
Press enter and you will see an automatically updating output message telling you the progress of the analysis. It should start like this:
Started analysis of ERR2935805.fastq
Approx 5% complete for ERR2935805.fastq
Approx 10% complete for ERR2935805.fastq
Approx 15% complete for ERR2935805.fastq
Approx 20% complete for ERR2935805.fastq
Approx 25% complete for ERR2935805.fastq
Approx 30% complete for ERR2935805.fastq
In total, it should take around ten minutes for FastQC to run on this fastq file (however, this depends on the size and number of files you give it). When the analysis completes, your prompt will return. So your screen will look something like this:
Approx 75% complete for ERR2935805.fastq
Approx 80% complete for ERR2935805.fastq
Approx 85% complete for ERR2935805.fastq
Approx 90% complete for ERR2935805.fastq
Approx 95% complete for ERR2935805.fastq
Analysis complete for ERR2935805.fastq
$
The FastQC program has created two new files within our
analysis/illumina_qc/ directory. We can see them by listing the contents of the illumina_qc folder
ls illumina_qc/
ERR2935805_fastqc.html ERR2935805_fastqc.zip
For each input FASTQ file, FastQC has created a .zip file and a .html file. The .zip file extension indicates that this is actually a compressed set of multiple output files. A summary report for our data is in the he .html file.
You need to transfer ERR2935805_fastqc.html from your AWS instance to your local computer to view it with a web browser.
To do this we will use the scp (secure copy protocol) command. You need to start a second terminal window that is not logged into the cloud instance and ensure you are in your cloudspan directory. This is important because it contains your pem file which will allow the scp command access to your AWS instance to copy the file.
Starting a new new terminal
Open your file manager and navigate to the
cloudspanfolder which should contain the login key file- Open your machine’s command line interface: Windows users: Right click anywhere inside the blank space of the file manager, then select Git Bash Here. Mac users: Open Terminal and type
cdfollowed by the absolute path that leads to yourcloudspanfolder. Press enter.- Check that you are in the right folder using
pwd
Now use scp to download the file.
The command will look something like:
scp -i login-key-instanceNNN.pem csuser@instanceNNN.cloud-span.aws.york.ac.uk:~/cs_course/analysis/qc/illumina_qc/ERR2935805_fastqc.html .
Remember to replace NNN with your instance number.
As the file is downloading you will see an output like:
scp -i login-key-instanceNNN.pem csuser@instanceNNN.cloud-span.aws.york.ac.uk:~/cs_course/analysis/qc/illumina_qc/ERR2935805_fastqc.html .
ERR2935805_fastqc.html 100% 591KB 1.8MB/s 00:00
Once the file has downloaded File Explorer (Windows) or Finder (Mac) to find the file and open it - it should open up in your browser.
Help!
If you had trouble downloading and viewing the file you can view it here: ERR2935805_fastqc.html
First we will look at the “Per base sequence quality” graph.
The x-axis displays the base position in the read, and the y-axis shows quality scores. In this example, the sample contains reads that are 202 bp long.
Each position has a box-and-whisker plot showing the distribution of quality scores for all reads at that position.
- The horizontal red line indicates the median quality score.
- The yellow box shows the 1st to 3rd quartile range (this means that 50% of reads have a quality score that falls within the range of the yellow box at that position).
- The whiskers show the absolute range, which covers the lowest (0th quartile) to highest (4th quartile) values.
The plot background is also color-coded to identify good (green), acceptable (yellow), and bad (red) quality scores.
In this sample, the quality values do not drop much lower than 32 at any position. This is a high quality score meaning the sequence is high quality. This means that we do not need to do any filtering. Lucky us!
We should also have a look at the “Adapter Content” graph which will show us where adapter sequences occur in the reads. Adapter sequences are short sequences that are added to the sample to aid during the preparation of the DNA library. They therefore don’t tell us anything biologically important and should be removed if they are present in high numbers. They might also be removed in the case of certain applications, such as ones when the base sequence needs to be particularly accurate.
This graph shows us that this sequencing file has a low percentage (~2-3%) of adapter sequences in the reads, which means we do not need to trim any adapter sequences either.
When sequencing is poor(er) Quality
While the sequencing in this example is high quality this will not always be the case.
Here is an example of a good quality FastQC output and a bad quality FastQC output. The programe cutadapt can be used to filter poor quality reads and trim poor quality bases. See Genomics - Trimming and Filtering to learn more about trimming and filtering poor quality reads.
Nanopore quality control
Next we will assess the quality of the Nanopore raw reads. These are found in the file located at ~/cs_course/data/nano_fastq/ERR3152367_sub5.fastq.
Let us again view the first complete read in one of the files from our dataset by using head to look at the first four lines.
Move to the folder containing the Nanopore data:
cd ~/cs_course/data/nano_fastq/
Use head to look at the first four line of the fastq file:
head -n 4 ERR3152367_sub5.fastq
@ERR3152367.34573250 d8c83b24-b46e-4f1a-836f-768f835acf68 length=320
GGTTGGTTATGTGCATGTTTTCAGTTACATATTGCATCTGTGGGAGCATATTCTTGTTTATGGGTTATGTGTTGGTGGTTGCATGTGGTGTGTTGTTGTGTTAACAAGTGTGGAACCTGTTCATTGGGTTATGAACAACGACACAAGTGTTGCGTGTTGAGCTAGTTAACGTGTGTGTTGTTATTCTTCTGAACCAGTTAACTTATTTGTTTTGTTGGGTGTGAAGCAGTGGGCGTGAAGGTGAGCGATGAAGCGGCGTTGTTCTGTTGCGTGTTTGATTGTGTTGTGTTGCGTGAAGAAGCGTCGTTGTTGGGTGGTTC
+
$$##$$###$#%###%##$%%$$###$#$$#$%##%&$$$$$$%#$$$$#$%#%$##$#$%#%$$#$$$%#$$#$%$$$$$#$%#$#$%$$$##$%%#&$#$#$$$$$%$$%$$%%$$#"$#$$$#&$$$$$#$$$$$######$#$#$$###$%###$$$$%$$&%$$$#$#$$%#%$##$##%#$&$$$$$#$$$%$$$##%#%$##$%%$$#$$$$%#%$###$$$####%$%%$$'$$%$$$$$%$#$$&$$%$#$##$%%$$%$$%%$%&'##$##%$#$$%$###$$$$$#$$$$#$&%##$$#$$%$$$%###
This read is 320 bp, longer than the Illumina reads we looked at earlier. The length of a raw read from Nanopore sequencing varies depends on the length of the length of the DNA strand being sequenced.
Line 4 shows us the quality score of this read.
$$##$$###$#%###%##$%%$$###$#$$#$%##%&$$$$$$%#$$$$#$%#%$##$#$%#%$$#$$$%#$$#$%$$$$$#$%#$#$%$$$##$%%#&$#$#$$$$$%$$%$$%%$$#"$#$$$#&$$$$$#$$$$$######$#$#$$###$%###$$$$%$$&%$$$#$#$$%#%$##$##%#$&$$$$$#$$$%$$$##%#%$##$%%$$#$$$$%#%$###$$$####%$%%$$'$$%$$$$$%$#$$&$$%$#$##$%%$$%$$%%$%&'##$##%$#$$%$###$$$$$#$$$$#$&%##$$#$$%$$$%###
Based on the PHRED quality scores (see above for a reminder) we can see that the quality score of the bases in this read are between 1-10, which is lower than the Illumina sequencing above.
Instead of using FastQC we will use a program called NanoPlot, which is installed on the instance, to create some plots for the whole sequencing file. NanoPlot is specially built for Nanopore sequences.
Other programs for Nanopore QC
Another popular program for QC of Nanopore reads is PycoQC.
It produces similar plots to NanoPlot but will also give you information about the overall Nanopore sequencing run. In order to generate these, PycoQC uses a
sequencing summaryfile generated by the Nanopore sequencer (e.g. MiniION or PromethION).We are using NanoPlot because the
sequencing summarythat PycoQC needs is not avaiable for this dataset. You can see example output files from PycoQC here: Guppy-2.1.3_basecall-1D_DNA_barcode.html.
We first need to navigate to the qc directory we made earlier cs_course/analysis/qc.
cd ~/cs_course/analysis/qc/
We are now going to run NanoPlot with the raw Nanopore sequencing file.
First we can look at the help documenation for NanoPlot to see what options are available.
NanoPlot --help
NanoPlot Help Documentation
We will use four flags when we run the NanoPlot command:
We also use --outdir to specify an output directory. We’re also going to use the flag --loglength to produce plots with a log scale and finally we’re going to use --threads to run the program on more than one thread to speed it up.
-
--fastqto specify the filetype and file to analyse. The raw Nanopore data is in the location/cs_workshop/data/nano_fastq/ERR3152367_sub5.fastqand we will use this full absolute path in the NanoPlot command. -
--outdirto specify the where the output files should be written. We are going to specifynano_qcso that NanoPlot will create a new directory in our current directory (qc) and write its output files to it. Note: with NanoPlot you don’t need to create this directory before running the command. -
--threadsspecifies how many threads to run the program on (more threads = more compute power = faster). We will specify 4 to indicate that four threads should be used. -
--loglengthspecifies that we want plots with a log scale.
NanoPlot --fastq ~/cs_course/data/nano_fastq/ERR3152367_sub5.fastq --outdir nano_qc --threads 4 --loglength
Now we have the command set up we can press enter and wait for NanoPlot to finish.
This will take a couple of minutes. You will know it is finished once your cursor has returned (i.e. you can type in the terminal again).
Once NanoPlot has finished we can have a look at the output.
First we need to navigate into the nano_qc directory NanoPlot created, then list the files.
cd nano_qc
ls
LengthvsQualityScatterPlot_dot.html LengthvsQualityScatterPlot_loglength_kde.png Non_weightedLogTransformed_HistogramReadlength.png
LengthvsQualityScatterPlot_dot.png NanoPlot_20221005_1630.log WeightedHistogramReadlength.html
LengthvsQualityScatterPlot_kde.html NanoPlot-report.html WeightedHistogramReadlength.png
LengthvsQualityScatterPlot_kde.png NanoStats.txt WeightedLogTransformed_HistogramReadlength.html
LengthvsQualityScatterPlot_loglength_dot.html Non_weightedHistogramReadlength.html WeightedLogTransformed_HistogramReadlength.png
LengthvsQualityScatterPlot_loglength_dot.png Non_weightedHistogramReadlength.png Yield_By_Length.html
LengthvsQualityScatterPlot_loglength_kde.html Non_weightedLogTransformed_HistogramReadlength.html Yield_By_Length.png
We can see that NanoPlot has generated a lot of different files.
Like before, we can’t view most of these files in our terminal as we can’t open images or HTML files. Instead we’ll download the core information to our own computer.
Luckily, the NanoPlot-report.html file contains all of the plots and information held in the other files so we only need to download that one onto our local computer using scp.
Use a terminal that is not logged into the cloud instance and ensure you are in your cloudspan directory. You may have one from earlier. If you do not, use the instructions above to start one.
Use scp to copy the file - the command will look something like:
scp -i login-key-instanceNNN.pem csuser@instanceNNN.cloud-span.aws.york.ac.uk:~/cs_course/analysis/qc/nano_qc/NanoPlot-report.html .
Remember to replace NNN with the instance number specific to you. As the file is downloading you will see an output like:
scp -i login-key-instanceNNN.pem csuser@instanceNNN.cloud-span.aws.york.ac.uk:~/cs_course/analysis/qc/nano_qc/NanoPlot-report.html .
NanoPlot-report.html 100% 3281KB 2.3MB/s 00:01
Once the file has downloaded File Explorer (Windows) or Finder (Mac) to find the file and open it - it should open up in your browser.
Help!
If you had trouble downloading and viewing the file you can view it here: NanoPlot-report.html
In the report we can view summary statistics followed by plots showing the distribution of read lengths and the read length vs average read quality.
Looking at the summary statistics table answer the following questions:
Exercise 1:
- How many sequences are in this file?
- How many bases are there in this entire file?
- What is the length of the longest read in the file and its associated mean quality score?
Solution
Quality Encodings Vary
Note that not all sequencing machines use the same encoding for quality. So
#might not always mean 3, a poor quality score.This means it’s essential that you know which sequencing platform was used to generate your data, so that you can tell your quality control program which encoding to use. If you choose the wrong encoding, you run the risk of throwing away good reads or (even worse) not throwing away bad reads! Nanopore quality encodings are no exception. You can read more about the differences with Nanopore sequencing here: EPI2ME - Quality Scores.
N50
The N50 length is a useful statistic when looking at sequences of varying length as it indicates that 50% of the total sequence is in reads (i.e. chunks) that are that size or larger.
For this FASTQ file 50% of the total bases are in reads that have a length of 5,373 bp or longer.
See the webpage What’s N50? for a good explanation. We will be coming back to this statistic in more detail when we get to the assembly step.
We can also look at some of the plots produced by NanoPlot.
One useful plot is the plot titled “Read lengths vs Average read quality plot using dots after log transformation of read lengths”.
This plot shows the average quality of the sequence against the read lengths. We can see that the majority of the sequences have a quality score at least 4, and low quality scores come from very short reads. This means that for this dataset we should remove those with a lower quality score in order to improve the overall quality of the raw sequences before assembling the metagenome.
Filtering Nanopore sequences by quality
We can use the program Seqkit (which contains many tools for FASTQ/A file manipulation) to filter our reads. We will be using the command seqkit seq to create a new file containing only the sequences with an average quality above a certain value.
After returning to our home directory, we can view the seqkit seq help documentation with the following command:
cd ~/cs_course/
seqkit seq -h
Seqkit seq help documentation
From this we can see that the flag -Q will “only print sequences with average quality qreater or equal than this limit (-1 for no limit) (default -1)”.
From the plot above we identified that many of the lower quality reads below 4 were shorter more here so we should set the minimum limit to 4.
seqkit seq -Q 4 data/nano_fastq/ERR3152367_sub5.fastq > data/nano_fastq/ERR3152367_sub5_filtered.fastq
In the command above we use redirection (>) to generate a new file data/nano_fastq/ERR3152367_sub5_filtered.fastq containing only the reads with an average quality of 4 or above.
We can now re-run NanoPlot on the filtered file to see how it has changed.
cd analysis/qc/
NanoPlot --fastq ~/cs_course/data/nano_fastq/ERR3152367_sub5_filtered.fastq --outdir nano_qc_filt --threads 4 --loglength
Once again, wait for the command to finish and then scp the NanoPlot-report.html to your local computer.
Help!
If you had trouble downloading the file you can view it here: NanoPlot-filtered-report.html
Compare the NanoPlot statistics of the Nanopore raw reads before filtering and after filtering and answer the questions below.
Exercise 2:
- How many reads have been removed by filtering?
- How many bases have been removed by filtering?
- What is the length of the new longest read and its associated average quality score?
Solution
Key Points
Quality encodings vary across sequencing platforms.
It is important to know the quality of our data to be able to make decisions in the subsequent steps.
Data cleaning is essential at the beginning of metagenomics workflows.
Due to differences in the sequencing technology Nanopore data must be handled differently.