Evaluating your raw sequencing data

Before you start the alignment and analysis processes, it us useful to perform some initial quality checks on your raw data. Here we will assume you have data from GSAF's Illumina HiSeq sequencer.

idev -m 60 -q development

Illumina sequence data format (FASTQ)

GSAF gives you paired end sequencing data in two matching fastq format files, contining reads for each end sequenced -- for example Sample_ABC_L005_R1.cat.fastq and Sample_ABC_L005_R2.cat.fastq. Each read end sequenced is representd by a 4-line entry in the fastq file.

A 4-line fastq file entry looks like this:

@HWI-ST1097:104:D13TNACXX:4:1101:1715:2142 1:N:0:CGATGT
GCGTTGGTGGCATAGTGGTGAGCATAGCTGCCTTCCAAGCAGTTATGGGAG
+
=<@BDDD=A;+2C9F<CB?;CGGA<<ACEE*1?C:D>DE=FC*0BAG?DB6

1. Line 1 is the read identifier, which describes the machine, flowcell, cluster, grid coordinate, end and barcode for the read. Except for the barcode information, read identifiers will be identical for corresponding entries in the R1 and R2 fastq files.
2. Line 2 is the sequence reported by the machine.
3. Line 3 is always '+' from GSAF (it can optionally include a sequence description)
4. Line 4 is a string of Ascii-encoded base quality scores, one character per base in the sequence. For each base, an integer quality score = -10 log(probabilty base is wrong) is calculated, then added to 33 to make a number in the Ascii printable character range.

See the Wikipedia FASTQ format page for more information.

Exercise: Examine the 2nd sequence in a FASTQ file

What is the 2nd sequence in the file /corral-repl/utexas/BioITeam/ngs_course/intro_to_mapping/data/SRR030257_1.fastq?

head /corral-repl/utexas/BioITeam/ngs_course/intro_to_mapping/data/SRR030257_1.fastq

Counting sequences

One of the first thing to check is that your fastq files are the same length, and that length is evenly divisible by four. The wc command (word count) using the -l switch to tell it to count lines, not words, is perfect for this:

wc -l $BI/ngs_course/intro_to_mapping/data/SRR030257_1.fastq

Exercise: Counting FASTQ file lines

How many sequences are in the FASTQ file above?

What if your fastq file has been compressed, for example by gzip? You can still count the lines, and you don't have to uncompress the file to do it:

gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz | wc -l

Here you use gunzip -c to write decompressed data to standard output (-c means "to console", and leaves the original .gz file untouched). You then pipe that output to wc -l to get the line count.

Exercise: Counting compressed FASTQ lines

How many sequences are in the compressed FASTQ file above?

echo $((2368720 / 4))

FASTQ Quality Assurance tools

The first order of business after receiving sequencing data should be to check your data quality. This often-overlooked step helps guide the manner in which you process the data, and can prevent many headaches.

FastQC

FastQC is a tool that produces a quality analysis report on FASTQ files.

Useful links:

• FastQC report for a good Illumina dataset
• FastQC report for a bad Illumina dataset
• Online documentation for each FastQC report

First and foremost, the FastQC "Summary" should generally be ignored. Its "grading scale" (green - good, yellow - warning, red - failed) incorporates assumptions for a particular kind of experiment, and is not applicable to most real-world data. Instead, look through the individual reports and evaluate them according to your experiment type.

The FastQC reports I find most useful are:

1. The Per base sequence quality report, which can help you decide if sequence trimming is needed before alignment.
2. The Sequence Duplication Levels report, which helps you evaluate library enrichment / complexity. But note that different experiment types are expected to have vastly different duplication profiles.
3. The Overrepresented Sequences report, which helps evaluate adapter contamination.

Running FastQC

FastQC is not currently available from the TACC module system, but the command-line version has been installed in the $BI/bin/FastQC directory (downloaded from the Babraham Bioinformatics web site; interactive GUI versions are also available for Windows and Macintosh).

FastQC creates a sub-directory for each analyzed FASTQ file, so we should copy the file we want to look at locally first. Here's how to run FastQC using the version we installed:

# setup
cds
mkdir fastqc_test
cd fastqc_test
cp $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz .

# running the program
$BI/bin/FastQC/fastqc Sample_Yeast_L005_R1.cat.fastq.gz

Exercise: FastQC results

What did FastQC create?

drwxrwxr-x 4 abattenh G-803889     4096 May 20 22:59 Sample_Yeast_L005_R1.cat_fastqc
-rw-rw-r-- 1 abattenh G-803889   198239 May 20 22:59 Sample_Yeast_L005_R1.cat_fastqc.zip
-rwxr-xr-x 1 abattenh G-803889 51065629 May 20 22:59 Sample_Yeast_L005_R1.cat.fastq.gz

Looking at FastQC output

You can't run a web browser directly from your "dumb terminal" command line environment. The FastQC results have to be placed where a web browser can access them. We put a copy at this URL:

http://loving.corral.tacc.utexas.edu/bioiteam/yeast_stuff/Sample_Yeast_L005_R1.cat_fastqc/fastqc_report.html

Exercise: Should we trim this data?

Based on this FastQC output, should we trim this data?

Samstat

The samstat program can also produce a quality report for FASTQ files, and it can also report on aligned sequences in a BAM file.

Again, this program is not available through the TACC module system but is available in our $BI/bin directory (which is on your $PATH because of our common profile). You should be able just to type samstat and see some documentation.

Running samstat on FASTQ files