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 even if you do), you may notice later that something looks fishy in the the output: for example, many of your reads are not mapping or the ends of many of your reads do not align. Both can give you clues about whether you need to process the reads to improve the quality of data that you are putting into your analysis.
Here we will assume you have data from GSAF's Illumina HiSeq sequencer.
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).
- Diagnose common issues in FASTQ read files that will negatively impact analysis.
- Trim adaptor sequences and low quality regions from the ends of reads to improve analysis.
Table of Contents
| Table of Contents |
|---|
| Warning | ||
|---|---|---|
When following along here, please start an idev session for running any example commands:
|
Illumina sequence data format (FASTQ)
GSAF gives you paired end sequencing data in two matching FASTQ format files, containing 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.
...
| Expand | ||
|---|---|---|
| ||
Executing the command above reports that the 2nd sequence has ID = |
Counting sequences
If you get an error from running a program, one of the first thing to check is that the length of your FASTQ files is evenly divisible by four and — if the program expects paired reads — that the R1 and R2 files have the same number of reads. The wc command (word count) using the -l switch to tell it to count lines, not words, is perfect for this:
...
| Expand | |||||
|---|---|---|---|---|---|
| |||||
The bash shell has a really strange syntax for arithmetic: it uses a double-parenthesis operator. Go figure.
|
FASTQ Evaluation 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 that could require you to redo an entire analysis after they rear their ugly heads.
FastQC
FastQC is a tool that produces a quality analysis report on FASTQ files.
...
| Expand | ||
|---|---|---|
| ||
|
Running FastQC
FastQC is available from the TACC module system on lonestar. Interactive GUI versions are also available for Windows and Macintosh and can be downloaded from the Babraham Bioinformatics web site.
...
| Expand | |||||
|---|---|---|---|---|---|
| |||||
The Sample_Yeast_L005_R1.cat.fastq.gz file is what we analyzed, so FastQC created the other two items. Sample_Yeast_L005_R1.cat_fastqc is a directory (the "d" in "drwxrwxr-x"), so use ls Sample_Yeast_L005_R1.cat_fastqc to see what's in it. Sample_Yeast_L005_R1.cat_fastqc.zip is just a Zipped (compressed) version of the whole directory. |
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. You should copy the results back to your local machine (via scp or a GUI secure ftp client) to open them in a web browser.
...
| Expand | ||
|---|---|---|
| ||
The Per base sequence quality report does not look good. The data should probably be trimmed (to a constant 40 or 50 bp) before alignment. |
Samstat
The samstat program can also produce a quality report for FASTQ files. (We also use it again later to report on aligned sequences in a BAM file).
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
| Code Block | ||
|---|---|---|
| ||
# setup cds mkdir samstat_test cd samstat_test cp $BI/gva_course/mapping/data/SRR030257_1.fastq . # run the program samstat SRR030257_1.fastq |
...
| Code Block | ||
|---|---|---|
| ||
http://loving.corral.tacc.utexas.edu/bioiteam/SRR030257_1.fastq.html |
FASTQ Processing Tools
Trimming low quality bases
Low quality base reads from the sequencer can cause an otherwise mappable sequence not to align. 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.
FASTX Toolkit
The FASTX-Toolkit provides a set of command line tools for manipulating fasta and fastq files. The available modules are described on their website. They include a fast fastx_trimmer utility for trimming fastq sequences (and quality score strings) before alignment.
...
| Expand | |||||
|---|---|---|---|---|---|
| |||||
Type fastx_ then tab to see their names
|
Adapter trimming
Data from RNA-seq or other library prep methods that resulted in very short fragments can cause problems with moderately long (100-250 base) reads since the 3' end of sequence can extend through to the 3' adapter at a variable position and even past the end of the fragment. This 3' adapter contamination can cause the "real" insert sequence not to align because the adapter sequence does not correspond to the bases at the 3' end of the reference genome sequence.
...
The GSAF website describes the flavaors of Illumina adapter and barcode sequence in more detail https://wikis.utexas.edu/display/GSAF/Illumina+-+all+flavors
Cutadapt
The cutadapt program is an excellent tool for removing adapter contamination. The program is not available through TACC's module system but we've installed a copy in our $BI/bin directory.
The most common application of cutadapt is to remove adapter contamination from small RNA library sequence data, so that's what we'll show here. Note that this step is increasingly needed for genomic sequencing of MiSeq data with 250 base reads.
Running cutadapt on small RNA library data
When you run cutadapt you give it the adapter sequence to trim, and this is different for R1 and R2 reads.
...
| Expand | |||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| |||||||||||||||||||||||||
Please refer to https://wikis.utexas.edu/display/GSAF/Illumina+-+all+flavors for Illumina library adapter layout. The top strand, 5' to 3', of a read sequence looks like this.
The -a argument to cutadapt is documented as the "sequence of adapter that was ligated to the 3' end". So we care about the <Read 2 primer> for R1 reads, and the <Read 1 primer> for R2 reads. The "contaminent" for adapter trimming will be the <Read 2 primer> for R1 reads. There is only one Read 2 primer:
The "contaminant" for adapter trimming will be the <Read 1 primer> for R2 reads. However, there are three different Read 1 primers, depending on library construction:
Since R2 reads are the reverse complement of R1 reads, the R2 adapter contaminent will be the RC of the Read 1 primer used. For ChIP-seq libraries where reads come from both DNA strands, the TruSeq Read 1 primer is always used.
For RNAseq libraries, we use the small RNA sequencing primer as the Read 1 primer.
|
Flexbar
Flexbar provides a flexible suite of commands for demultiplexing barcoded reads and removing adapter sequences or low quality regions from the ends of reads.
...
Note that flexbar only searches for the sequences given (with options to allow for a given number of mismatches) NOT the reverse complement of those sequences therefore you must provide them yourself.
Trimmomatic
Trimmomatic offers similar options to Flexbar with the potential benefit that many illumina adaptor sequences are already "built-in". It is available here.