Quality control and filtering of the raw sequence files

Authors
Affiliation

Jennifer Lu

EMBL-EBI

Varsha Kale

EMBL-EBI

Published

August 30, 2023

Prerequisites

These instructions are for the course VM. To run externally, please refer to the section at the end.

For this tutorial, you’ll need to move into the working directory and start a Docker container. Set the variable DATADIR as instructed.

cd /home/training/quality
chmod -R 777 /home/training/quality
export DATADIR=/home/training/quality
xhost +

You will get the message “access control disabled, clients can connect from any host” Now start the Docker container:

docker run --rm -it  -e DISPLAY=$DISPLAY  -v $DATADIR:/opt/data -v /tmp/.X11-unix:/tmp/.X11-unix:rw -e DISPLAY=unix$DISPLAY microbiomeinformatics/biata-qc-assembly:v2021

Quality control and filtering of the raw sequence files

Learning Objectives

In the following exercises, you’ll learn how to check the quality of short read sequences, identify adaptor sequences, remove adapters and low-quality sequences, and construct a reference database for host decontamination.

Here you should see the contents of the working directory.

These are the files we’ll use for the practical. Move into the folder:

ls /opt/data
cd /opt/data
Generate a directory of the FastQC results 
mkdir fastqc_results
fastqc oral_human_example_1_splitaa.fastq.gz
fastqc oral_human_example_2_splitaa.fastq.gz
mv *.zip fastqc_results
mv *.html fastqc_results
chown 1001 fastqc_results/*.html
Now on your computer, select the folder icon.

Navigate to Home → quality → fastqc_results

Right-click on file oral_human_example_1_splitaa_fastqc.html, select ‘open with other application’, and open with Firefox.

Screenshot of fast qc

Spend some time looking at the ‘Per base sequence quality.’

For each position, a BoxWhisker-type plot is drawn:

  • The central red line is the median value.
  • The yellow box represents the inter-quartile range (25-75%).
  • The upper and lower whiskers represent the 10% and 90% points.
  • The blue line represents the mean quality.

The y-axis on the graph shows the quality scores. The higher the score, the better the base call. The background of the graph divides the y-axis into very good quality calls (green), calls of reasonable quality (orange), and calls of poor quality (red). The quality of calls on most platforms will degrade as the run progresses, so it’s common to see base calls falling into the orange area towards the end of a read.

Question

What does this tell you about your sequence data? When do the errors start?

In the pre-processed files, we see two warnings, as shown on the left side of the report. Navigate to the “Per bases sequence content.”

Screenshot of fast qc
Question

At around 15-19 nucleotides, the DNA composition becomes very even; however, at the 5’ end of the sequence, there are distinct differences. Why do you think that is?

Step

Open up the FastQC report corresponding to the reversed reads.

Question

Are there any significant differences between the forward and reverse files?

For more information on the FastQC report, please consult the ‘Documentation’ available from this site: FastQC Documentation

We are currently only looking at two files, but often we want to look at many files. The tool multiqc aggregates the FastQC results across many samples and creates a single report for easy comparison. Here we will demonstrate the use of this tool.

Run 
cd /opt/data
mkdir multiqc_results
multiqc fastqc_results -o multiqc_results
chown 1001 multiqc_results/*.html

In this case, we provide the folder containing the FastQC results to multiqc, and the -o allows us to set the output directory for this summarized report.

Now on your computer, select the folder icon.

Navigate to Home → quality → multiqc_results

Right-click on file multiqc_report.html, select ‘open with other application’, and open with Firefox.

Screenshot of multiQC

Scroll down through the report. The sequence quality histograms show the above results from each file as two separate lines. The ‘Status Checks’ show a matrix of which samples passed check and which ones have problems.

Question

What fraction of reads are duplicates?

So far we have looked at the raw files and assessed their content, but we have not done anything about removing duplicates, sequences with low quality scores, or removal of the adaptors. So, let’s start this process. The first step in the process is to make a database relevant for decontaminating the sample. It is always good to routinely screen for human DNA (which may come from the host and/or staff performing the experiment). However, if the sample is from a mouse, you would want to download the mouse genome.

In the following exercise, we are going to use two “genomes” already downloaded for you in the decontamination folder. To make this tutorial quicker and smaller in terms of file sizes, we are going to use PhiX (a common spike in) and just chromosome 10 from human.

Run 
cd /opt/data/decontamination

For the next step, we need one file, so we want to merge the two different fasta files. This is simply done using the command-line tool cat.

Run 
cat phix.fasta GRCh38_chr10.fasta > GRCh38_phix.fasta

Now we need to build a bowtie index for them:

Run 
bowtie2-build GRCh38_phix.fasta GRCh38_phix.index

It is possible to automatically download a pre-indexed human genome in Bowtie2 format using the following command (but DO NOT do this now, as this will take a while to download):

kneaddata_database –download human_genome bowtie2

Now we are going to use the GRCh38_phix database and clean up our raw sequences. kneaddata is a helpful wrapper script for a number of pre-processing tools, including Bowtie2 to screen out contaminant sequences, and Trimmomatic to exclude low-quality sequences. We also have written wrapper scripts to run these tools (see below), but using kneaddata allows for more flexibility in options.

Run 
cd /opt/data
mkdir clean

We now need to uncompress the fastq files.

Uncompress 
gunzip -c oral_human_example_2_splitaa.fastq.gz > oral_human_example_2_splitaa.fastq
gunzip -c oral_human_example_1_splitaa.fastq.gz > oral_human_example_1_splitaa.fastq

kneaddata --remove-intermediate-output -t 2 --input oral_human_example_1_splitaa.fastq --input oral_human_example_2_splitaa.fastq --output /opt/data/clean --reference-db /opt/data/decontamination/GRCh38_phix.index --bowtie2-options "--very-sensitive --dovetail" --trimmomatic /opt/data/Trimmomatic-0.39/ --bypass-trf --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50"

The options above are:

  • --input, Input FASTQ file. This option is given twice as we have paired-end data.
  • --output, Output directory.
  • --reference-db, Path to Bowtie2 database for decontamination.
  • -t, Number of threads to use (2 in this case).
  • --trimmomatic-options, Options for Trimmomatic to use, in quotations.
  • --bowtie2-options, Options for Bowtie2 to use, in quotations.
  • --remove-intermediate-output, Intermediate files, including large FASTQs, will be removed.

Kneaddata generates multiple outputs in the “clean” directory, containing different four different files for each read.

Run FastQC

Using what you have learned previously, generate a FastQC report for each of the oral_human_example_1_splitaa_kneaddata_paired files. Do this within the clean directory.

cd /opt/data/clean
mkdir fastqc_final
<you construct the commands>
mv /opt/data/clean/*.zip /opt/data/clean/fastqc_final
mv /opt/data/clean/*.html /opt/data/clean/fastqc_final
chown 1001 /opt/data/clean/fastqc_final/*.html
Run multiQC

Also generate a multiQC report. Send the output to the folder multiqc_final and look at the sequence quality histograms.

cd /opt/data/clean/
mkdir multiqc_final
<you construct the command>
chown 1001 /opt/data/clean/multiqc_final/*.html
Check report

View the MultiQC report as before using your browser.

Screenshot of multiQC

Scroll down through the report. The sequence quality histograms show the above results from each file as two separate lines. The ‘Status Checks’ show a matrix of which samples passed check and which ones have problems.

Question

Open the previous MultiQC report and see if they have improved?

Did sequences at the 5’ end become uniform? Why might that be? Is there anything that suggests that adaptor sequences were found?

To generate a summary file of how the sequences were categorized by Kneaddata, run the following command.

Run 
cd /opt/data/clean
kneaddata_read_count_table --input /opt/data/clean --output kneaddata_read_counts.txt
cat kneaddata_read_counts.txt
Question

What fraction of reads have been deemed to be contaminating?

The reads have now been decontaminated and can be uploaded to ENA, one of the INSDC members. It is beyond the scope of this course to include a tutorial on how to submit to ENA, but there is additional information available on how to do this in this Online Training guide provided by EMBL-EBI

Assembly PhiX decontamination

Learning Objectives

In the following exercises, you will generate a PhiX blast database and run a blast search with a subset of assembled freshwater sediment metagenomic reads to identify contamination.

PhiX, used in the previous section of this practical, is a small bacteriophage genome typically used as a calibration control in sequencing runs. Most library preparations will use PhiX at low concentrations; however, it can still appear in the sequencing run. If not filtered out, PhiX can form small spurious contigs that could be incorrectly classified as diversity.

Generate the PhiX reference blast database: 
cd /opt/data/decontamination
makeblastdb -in phix.fasta -input_type fasta -dbtype nucl -parse_seqids -out phix_blastDB

Prepare the freshwater sediment example assembly file and search against the new blast database. This assembly file contains only a subset of the contigs for the purpose of this practical.

Run 
cd /opt/data
gunzip -c freshwater_sediment_contigs.fa.gz > freshwater_sediment_contigs.fa
blastn -query freshwater_sediment_contigs.fa -db decontamination/phix_blastDB -task megablast -word_size 28 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 0.0001 -min_raw_gapped_score 100 -penalty -5 -soft_masking true -window_size 100 -outfmt 6 -out freshwater_blast_out.txt

The blast options are:

  • -query, Input assembly fasta file.
  • -out, Output file
  • -db, Path to blast database.
  • -task, Search type -“megablast”, for very similar sequences (e.g, sequencing errors)
  • -word_size, Length of initial exact match
Add headers to the blast output and look at the contents of the final output file: 
cat blast_outfmt6.txt freshwater_blast_out.txt > freshwater_blast_out_headers.txt
cat freshwater_blast_out_headers.txt
Question

Are the hits significant?

Question

What are the lengths of the matching contigs? We would typically filter metagenomic contigs at a length of 500bp. Would any PhiX contamination remain after this filter?

Now that PhiX contamination was identified, it is important to remove these contigs from the assembly file before further analysis or upload to public archives.

Using Negative Controls

Learning Objectives

This exercise will look at the analysis of negative controls. You will assess the microbial diversity between a negative control and a skin sample.

The images below show the taxonomic classification of two samples: a reagent negative control and a skin metagenomic sample. The skin sample is taken from the antecubital fossa - the elbow crease, which is moist and site of high microbial diversity. The classification was performed with kraken2. Kraken2 takes a while to run, so we have done this for you and plotted the results. An example of the command used to do this. DO NOT run this now:

kraken2 --db standard_db --threshold 0.10 --threads 8 --use-names --fastq-input --report out.report --gzip-compressed in_1.fastq.gz in_2.fastq.gz See the kraken2 manual for more information

See Pavian manual for the plots.

The following image shows the microbial abundance in the negative control: Kraken negative control

The following image shows the microbial abundance in the skin sample: Kraken skin sample

Step

Look for similarities and differences at both the phylum and genus level - labelled as ‘P’ and ‘G’ on the bottom axis.

Question

Is there any overlap between the negative control and skin sample phylum? Can we map the negative control directly to the skin sample to remove all contaminants? If not, why?

Question

Are there any genera in the negative control which aren’t present in the skin sample? If you do a google search of this genus, where are they commonly found? With this information, where could this bacteria in the negative control have originated from?

If you have finished the practical you can try this step for more practice assessing and trimming datasets, there is another set of raw reads called “skin_example_aa” from the skin metagenome available. These will require a fastqc or multiqc report, followed by trimming and mapping to the reference database with kneaddata. Using what you have learned previously, construct the relevant commands. Remember to check the quality before and after trimming.

Tip

Consider other trimmomatic options from the manual e.g. “ILLUMINACLIP”, where /opt/data/NexteraPE-PE is a file of adapters.

Navigate to skin folder and run quality control. 
cd /opt/data/skin
<construct the required commands>

Remember you will need to run the following command to view any html files in the VM browsers:

chown 1001 foldername/*.html

Running the practical externally

We need to first fetch the practical datasets.

    wget http://ftp.ebi.ac.uk/pub/databases/metagenomics/mgnify_courses/ebi_2020/quality.tar.gz
    or
    rsync -av --partial --progress rsync://ftp.ebi.ac.uk/pub/databases/metagenomics/mgnify_courses/ebi_2020/quality.tar.gz .

Once downloaded, extract the files from the tarball:

    tar -xzvf quality.tar.gz

We also need the trimmomatic binary

    cd quality
    wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
    unzip Trimmomatic-0.39.zip
    cd ..

Now pull the docker container and set the above quality directory as DATADIR.

    docker pull microbiomeinformatics/biata-qc-assembly:v2021
    export DATADIR={path to quality directory}
    xhost +

You will see the message “access control disabled, clients can connect from any host”

    docker run --rm -it  -e DISPLAY=$DISPLAY  -v $DATADIR:/opt/data -v /tmp/.X11-unix:/tmp/.X11-unix:rw -e DISPLAY=unix$DISPLAY microbiomeinformatics/biata-qc-assembly:v2021

The container has the following tools installed: - kneaddata - fastqc - multiqc - blast - bowtie-2

You can now continue this practical from the section “Quality control and filtering of the raw sequence files”

Reuse

Apache 2.0