2. Quality control & removal of optical duplicates

1. How to use the exercise documents

The first thing to notice is that there are links to the different sections in the exercise in the left hand menu of each exercise together with a green box containng a short introduction to the exercise or tool you will be using

You will also find links to external web sites describing most of the tools you will be using throughout the exercises

You can copy and paste commands directly from this HTML and into your terminal if you don’t want to write the complete commands yourself

The commands are written in boxes like this:

any command

File paths to data are also written like the commands, for example: path/to/files

The tilde (~) symbol means home directory, therefore ~/practical/ is the same as /home/practical

We try to write tool in bold letters and filnames in boxes like this filename, but we are only human…

The tasks you will be performing will be annonced like this:

I] Read this line which is printed in blue

You will be given questions troughout the execises. The questions will be given like this:

❓ Can you read this?

The solution to the exercises is found below each question like this:

Hints, notes and useful information are (often) put in blue boxes like below:

Optional exercises for those that work fast or want to continue outside the course hours are put in yellow boxes like below.

2. Make sure all your data is there

I] In your home directory, create a directory named practical:

mkdir practical

II] Move to this directory:

cd practical

III] Copy the data for this exercise from the shared disk:

cp -r /net/software/practical/2 .

IV] By default you don’t have permission to write to the directory you copied from the shared disk, so you need to change the permissions:

chmod -R a+rwx 2

# chmod, change permission
# -R, recursive (to all subfolders and files in the current directory)
# a+rwx, give permission to read, write and execute
# 2, the directory you are changing the permission on

V] Move to the directory named rawdata (cd rawdata) and check that you are in the directory named ~/practical/2/rawdata:

pwd

VI] Check what data is in this directory:

l

VII] Open the md5_data.txt file and look at the content in a text editor. It should have one checksum per file.

VIII] Check that the sum is correct for your files by typing in the terminal:

md5sum -c md5_data.txt

If the files are complete, the test should return ‘OK’ for all files

Optional: You can try to generate MD5 checksum on a file and output it to a file:

md5sum filename > md5_test.txt

3. Get to know the FASTQ file format

I] In the directory where the compressed FASTQ files are located, type:

file sample1_R1.fastq.gz

❓ What type of file is this?

II] Try to run the following command:

zcat sample1_R1.fastq.gz | head

❓ Each sequence read has four lines. What information does these lines contain?

III] You can preview more of the file, eg. the 100 first lines by:

zcat sample1_R1.fastq.gz | head -n 100

It may be useful to know how many sequences a FASTQ file contain. Also it is useful to convert FASTQ files to FASTA format.

IV] Count numbers of sequence reads in the compressed file:

zcat sample1_R1.fastq.gz | paste - - - - | wc -l

❓ Are there equal number of sequence reads in both the paired FASTQ files (and how many reads)?

4. Simple FASTQ/FASTA manipulations using SeqKit

I] In the directory ~/practical/2/, create a directory named seqkit and move to the directory

II] Make a symbolic links to the raw data:

ln -s ~/practical/2/rawdata/sample1_R* .

III] Use Seqkit to produce simple statistics of your files with the subcommand stats:

seqkit stats *.fastq.gz

IV] Try to use the -a flag

❓ What is the main difference between the outputs?

❓ What is the average length of the sequence reads in sample1_R1.fastq.gz?

❓ How many percent of the sequence reads in sample1_R1.fastq.gz have an average above Q30?

V] View the top two sequence reads in the forward FASTQ file using the subcommand head:

seqkit head -n 2 sample1_R1.fastq.gz

VI] View the top two sequence reads in the forward FASTQ file using the subcommand range:

seqkit range -r 1:2 sample1_R1.fastq.gz

VII] View the last two sequence reads in the forward FASTQ file using the subcommand range:

seqkit range -r -2:-1 sample1_R1.fastq.gz

VIII] You can combine several subcommands, eg. to produce simple statistics for the first 100 sequence reads:

cat sample1_R1.fastq.gz | seqkit range -r 1:100 | seqkit stats -a

IX] Convert fastq.gz to FASTA using the subcommand fq2fa:

seqkit fq2fa sample1_R1.fastq.gz -o sample1_R1.fasta

gunzip -c sample1_R1.fastq.gz | awk '{if(NR%4==1) {printf(">%s\n",substr($0,2));} else if(NR%4==2) print;}' > sample1_R1.fasta

X] Have a look at the FASTA file:

head sample1_R1.fasta

❓ What is the main difference between a FASTQ and a FASTA file?

XI] Convert the reverse FASTQ to FASTA:

seqkit fq2fa sample1_R2.fastq.gz -o sample1_R2.fasta

XII] Concatenate the two FASTA files (head to tail) using cat and write out the new file:

cat sample1_R1.fasta sample1_R2.fasta > sample1_cat.fasta

It is possible to use Seqkit in order to subsample sequences by number or proportion.

XIII] Subsample 100 read pairs from the FASTQ files using the subcommand sample:

seqkit sample -n 100  -s 0 sample1_R1.fastq.gz -o sub_sample1_R1.fastq.gz
seqkit sample -n 100  -s 0 sample1_R2.fastq.gz -o sub_sample1_R2.fastq.gz

# -n, sample by number (result may not exactly match)
# -s, random seed

❓ How many sequence reads are there in the reduced file:

XIV] Subsample 10% of the reads from the R1 FASTQ file using the subcommand sample. Hint: the flag -p = sample by proportion.

❓ Count the number of sequence in one of the reduced files. How many sequence reads were there?

Optional exercises

If you remember from the lecture, the FASTQ header contains information from the sequencing machine, among other the position on the flow cell:

If you check the FASTQ headers in your sub_sample1_R1.fastq.gz you will notice that the sequence reads were produced on different tiles on the flow cell.

---

XV] Use the subcommand grep to extract all sequence reads from tile 1101. Hint: read the uncompressed file into the memory using zcat:

❓How many sequence reads were there from tile 1101 in sub_sample1_R1.fastq.gz?

Solution - Click to expand

5 reads. One way to find out is:

zcat sub_sample1_R1.fastq.gz | seqkit grep -n -r -p 1101 | seqkit stats

XVI] You can use the subcommand grep to extract all sequence reads that has a certain sequence pattern. Extract all reads containing the pattern AACTCAT from sub_sample1_R1.fastq.gz. Hint: read the uncompressed file into the memory using zcat:

❓How many sequence reads contained the pattern AACTCAT in sub_sample1_R1.fastq.gz?

Solution - Click to expand

3 reads

zcat sub_sample1_R1.fastq.gz | seqkit grep -s -i -p AACTCAT | seqkit stats

5. Quality control using FastQC

I] In the directory ~/practical/2/, create a directory named fastqc and move to the directory

II] Make a symbolic links to the raw data:

ln -s ~/practical/2/rawdata/sample1_R* .

III] Start FastQC by typing fastqc in a terminal. After a few second the graphical FastQC interface will appear.

VI] In the FastQC GUI, click on “File” and select the two compressed FASTQ files sample1_R1.fastq.gz and sample1_R2.fastq.gz in the ~/practical/2/rawdata directory.

❓ What does green tick, the orange triangle and the red cross indicate?

❓ How many sequence reads does each of your dataset consist of?

❓ What is the sequence length range?

❓ What is the GC content of the metagenomic sequence data?

❓ What quality format are your sequence reads in?

❓ What does a quality score of 30 indicate about the accuracy?

❓ Which of the sequence files would you say have the best base calling quality?

6. Remove optical duplicates from NGS data

I] In the directory ~/practical/2/, create a directory named clumpify and move to the directory

II] Make a symbolic links to the raw data:

ln -s ~/practical/2/rawdata/sample1_R* .

III] Remove optical duplicates and tile-edge duplicates. NB! YOU HAVE TO SET DISTANCE (replace the X in “dist=X”).

clumpify.sh in1=sample1_R1.fastq.gz in2=sample1_R2.fastq.gz out1=sample1_clumped_R1.fastq.gz out2=sample1_clumped_R2.fastq.gz dedupe optical dist=X -Xmx8g
		
# dedupe, Remove duplicates - combined with optical dist = Remove optical duplicates
# optical dist=X, duplicates within X pixels of each other
# -Xmx8g, set memory to 8 GB RAM

OUTPUT: 
sample1_clumped_R1.fastq.gz
sample1_clumped_R2.fastq.gz

IV] When Clumpify finishes, look at the number of reads and file size of the original FASTQ files and the “clumpified” files.

❓ What does the flag in1 do?

❓ What does the flag out1= do?

❓ How many duplicates were identified?

❓ How many percent smaller is the size of the “clumpified” file of R1 relative to the original file?

As you can see, even only a small number of reads were removed as being optical duplicates, the file compression is greatly improved by Clumpify. Gzip compression becomes much more efficient by reordering reads so that reads with similar sequence are nearby.