Overview

Before you start the alignment and analysis processes, it can be useful to perform some initial quality checks on your raw data. If you don't do this (or don't do this sufficiently), you may notice at the end of your analysis some things still are not clear: for example, maybe a large portion of reads do not map to your reference or maybe the reads map well, except the ends do not align at all. Both of these results can give you clues about how you need to process the reads to improve the quality of data that you are putting into your analysis.

Learning Objectives

This tutorial covers the commands necessary to use several common programs for evaluating read files in FASTQ format and for processing them (if necessary).

1. Use basic linux commands to determine read count numbers and pull out specific reads.
2. Diagnose common issues in FASTQ read files that will negatively impact analysis.
3. Trim adaptor sequences and low quality regions from the ends of reads to improve analysis.

FASTQ data format

A common question is 'after you submit something for sequencing what do you get back?' The answer is FASTQ files.

...

See the Wikipedia FASTQ format page for more information.

Determine 2nd sequence in a FASTQ file

What the 2nd sequence in the file $BI/gva_course/mapping/data/SRR030257_1.fastq is?

...

1. The -n option can be used to control how many lines of a file are printed:

   Code Block
   language bash title Show just the first 2 reads collapse true
   head -n 8 $BI/gva_course/mapping/data/SRR030257_1.fastq

2. The output of the head command can be piped to the tail command to isolate specific groups of lines:

   Code Block
   language bash title Show just the 2nd read collapse true
   head -n 8 $BI/gva_course/mapping/data/SRR030257_1.fastq | tail -n 4

3. The grep command can be used to look for lines that contain only ACTG or N:

   Code Block
   language bash title head | tail | grep to isolate the 2nd read collapse true
   head -n 8 $BI/gva_course/mapping/data/SRR030257_1.fastq | tail -n 4 | grep "^[ACTGN]*$"

   Code Block
   language bash title head to show the first 10 lines, grep to only print the 3 sequence lines collapse true
   head $BI/gva_course/mapping/data/SRR030257_1.fastq | grep "^[ACTGN]*$"

   Code Block
   language bash title grep to only print the first 3 lines containing sequence collapse true
   grep -m 2 "^[ACTGN]*$" $BI/gva_course/mapping/data/SRR030257_1.fastq

   ^ by increasing the "-m" value we can now quickly get a block of sequence of any size of our choosing. This is the first truly useful command on this page. With a block of sequence, you can start to see things like:

   1. if the first/last bases are always the same

   2. if the reads are the same length

   3. if a single sequence shows up a huge number of times

4. Tip
   This is our first example that there are often many different ways to do the same thing in NGS analysis. While some may be faster, or more efficient, the same answers are still achieved. Don't let the pursuit of perfection keep you from getting your answers in whatever way you can justify.

Counting sequences

Often, the first thing you (or your boss) want to know about your sequencing run is simply, "how many reads are there?". For the $BI/gva_course/mapping/data/SRR030257_1.fastq file, the answer is 3,800,180. How can we figure that out?

The grep (or Global Regular Expression Print) command can be used to determine the number of lines which match some criteria as shown above. Above we used it to search for:

1. anything from the group of ACTGN with the [] marking them as a group
2. matching any number of times *
3. from the beginning of the line ^
4. to the end of the line $

Here, since we are only interested in the number of reads that we have, we can make use of knowing the 3rd line in the fastq file is a + and a + only, and grep's -c option to simply report the number of reads in a file.

grep -c "^+$" $BI/gva_course/mapping/data/SRR030257_1.fastq

...

grep -c "^[ACTGN]*$" $BI/gva_course/mapping/data/SRR030257_1.fastq

wc -l $BI/gva_course/mapping/data/SRR030257_1.fastq 

echo $((15200720 / 4))

Counting with compressed files

Thus far we have looked at a FASTQ file. Because FASTQ files contain millions-billions of lines and billions+ characters, they are often stored in a compressed format. Specifically they are typically stored in a 'gzipped' format to save storage space and typically will end with ".gz" so you can identify them. While the files are easily changed from compressed to noncompressed and back again (and you will do some of this throughout the course and plenty more in your own work), the bigger the file, the longer such actions will take.

...

While checking the number of reads a file has can solve some of the most basic problems, it doesn't really provide any direct evidence as to the quality of the sequencing data. To get this type of information before starting meaningful analysis other programs must be used.

Evaluating FASTQ files with FastQC

FastQC overview

Once you move past the most basic questions about your data, you need to move onto more substantive assessments. As discussed above, this often-overlooked step helps guide the manner in which you process the data, and can prevent many headaches that could require you to redo an entire analysis after they rear their ugly heads.

...

Running FastQC

You may recall from today's first tutorial, that we used the conda system to install fastqc in preparation for this tutorial. If you did not complete that, please go back and do so now, and don't hesitate to ask a question if you are having difficulties. Interactive GUI versions are also available for Windows and Macintosh and can be downloaded from the Babraham Bioinformatics web site. We don't want to clutter up our work space so copy the SRR030257_1.fastq file to a new directory named GVA_fastqc_tutorial on scratch, and use fastqc's help option after the module it is loaded installed to figure out how to run the program. Once the program is completed use scp to copy the important file back to your local machine (The bold words are key words that may give you a hint of what steps to take next)

...

-rwxr-xr-x 1 ded G-802740 498588268 May 23 12:06 SRR030257_1.fastq
-rw-r--r-- 1 ded G-802740    291714 May 23 12:07 SRR030257_1_fastqc.html
-rw-r--r-- 1 ded G-802740    455677 May 23 12:07 SRR030257_1_fastqc.zip

Looking at FastQC output

As discussed in the introduction tutorial, you can't run a web browser directly from your command line environment. You should copy the results back to your local machine (via scp) so you can open them in a web browser.

...

# on tacc terminal
pwd
 
# on new terminal of local computer
scp <username>@ls5<username>@stampede2.tacc.utexas.edu:<pwd_results_from_other_window>/SRR030257_1_fastqc.html ~/Desktop
 
# open the newly transferred file from from the desktop and see how the data looks

...

FASTQ Processing Tools

Cutadapt

There are a number of open source tools that can trim off 3' bases and produce a FASTQ file of the trimmed reads to use as input to the alignment program. Cutadapt provides a simple command line tool for manipulating fasta and fastq files. The program description on their website provides good details of all the capabilities and examples for some common tasks. Using what you learned about installing fastqc via conda, see if you can figure out how to install cutadapt.

...

To run the program, you simply type 'cutadapt' followed by whatever options you want, and then the name of the fastq files without any option in front of it. Use the -h option to see all the different things you can do and what options you think might be useful for removing a string of "A"s at the 3` end of the read, and shortening the read to 34 bp

Adapter trimming

cutadapt can be used to trim specific sequences. Typically this is done to remove adapter sequences introduced during library prep that are not related to your sample. Based on our fastqc analysis, the sequence AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA is significantly overrepresented in our data. How would you use cutadapt to remove those sequences from the fastq file?

...

Trimming low quality bases

Low quality base reads from the sequencer can cause an otherwise mappable sequence not to align. See if you can come up with a cutadapt command to trim the reads down to 34bp.

...

cutadapt -o SRR030257_1.Adepleted_and_trimmed.fastq -l 34 -a AAAAAAAAAAAAAAAAAAAA -m 16 SRR030257_1.fastq

Bonus Exercise: compressing the trimmed files 

As mentioned above, fastq files are often stored as gzipped compressed files as they are smaller, easier to transfer, and many programs allow for their use directly.

...

This is another example of different solutions giving the same final product, and how careful reading of documentation may improve your work. NGS data analysis is a results driven process, if the result is correct, and you know how you got the result it is correct as long as it is reproducible.

A note on versions 

cutadapt --version
conda list
conda list cutadapt

Optional Exercise: Improve the quality of R2 the same way you did for R1.

Unfortunately we don't have time during the class to do this, but as a potential exercise in your free time, you could improve R2 the same way you did R1 and use the improved fastq files in the subsequent read mapping and variant calling tutorials to see the difference it can make in overall results.

...