Pre-processing FASTQ sequences

Before you start the alignment and analysis processes, it us useful to perform some initial quality checks on your raw data. You may also need to pre-process the sequences to trim them or remove adapters. Here we will assume you have paired-end data from GSAF's Illumina HiSeq sequencer.

Setup

Special stampede login

Before we start, log into stampede like you did yesterday, but use this special hostname:

login8.stampede.tacc.utexas.edu

Remember how we emphasized yesterday that you should not perform significant computation on login nodes? Well, there are a  few exceptions, and login8.stampede.tacc.utexas.edu is one of them. Is it a dedicated login node owned by an organization our lab belongs to, so we have given you access to it for the duration of this course. This will let us do a few things at the command line that would normally set off alarm bells from the TACC folks if we all did them on a standard login node.

Data staging

Recall that yesterday we copied some fastq data to our "permanent archive area".

mkdir -p $WORK/archive/original/2014_05.yeast
cd $WORK/archive/original/2014_05.yeast

wget http://web.corral.tacc.utexas.edu/BioITeam/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz

wget http://web.corral.tacc.utexas.edu/BioITeam/yeast_stuff/Sample_Yeast_L005_R2.cat.fastq.gz

Let's now set ourselves up to process this data in $SCRATCH, using some of best practices for organizing our workflow.

# Create a $SCRATCH area to work on data for this course,
# with a sub-directory for pre-processing raw fastq files
mkdir -p $SCRATCH/core_ngs/fastq_prep

# Make a symbolic link to the original yeast data in $WORK
cd $SCRATCH/core_ngs/fastq_prep
ln -s $WORK/archive/original/2014_05.yeast fq

Illumina sequence data format (FASTQ)

GSAF gives you paired end sequencing data in two matching fastq format files, containing reads for each end sequenced.

# make sure you're in the $SCRATCH/core_ngs/fastq_prep directory
ls fq  

Sample_Yeast_L005_R1.cat.fastq.gz Sample_Yeast_L005_R2.cat.fastq.gz

4-line FASTQ format

Each read end sequenced is represented by a 4-line entry in the fastq file that looks like this:

@HWI-ST1097:127:C0W5VACXX:5:1101:4820:2124 1:N:0:CTCAGA
TCTCTAGTTTCGATAGATTGCTGATTTGTTTCCTGGTCATTAGTCTCCGTATTTATATTATTATCCTGAGCATCATTGATGGCTGCAGGAGGAGCATTCTC
+
CCCFFFFDHHHHHGGGHIJIJJIHHJJJHHIGHHIJFHGICA91CGIGB?9EF9DDBFCGIGGIIIID>DCHGCEDH@C@DC?3AB?@B;AB??;=A>3;;  

Line 1 is the unique read name. The format is as follows, using the read name above:

machine_id:lane:flowcell_grid_coordinates  end_number:failed_qc:0:barcode

@HWI-ST1097:127:C0W5VACXX:5:1101:4820:2124 1:N:0:CTCAGA

• The line as a whole will be unique for this read fragment. However, the corresponding R1 and R2 reads will have identical machine_id:lane:flowcell_grid_coordinates information. This common part of the name ties the two read ends together.
• Most sequencing facilities will not give you qc-failed reads (failed_qc = Y) unless you ask for them.

Line 2 is the sequence reported by the machine, starting with the first base of the insert (the 5' adapter has been removed by the sequencing facility). These are ACGT or N uppercase characters.

Line 3 is always '+' from GSAF (it can optionally include a sequence description)

Line 4 is a string of Ascii-encoded base quality scores, one character per base in the sequence.

• For each base, an integer Phred-type quality score is calculated as integer score = -10 log(probabilty base is wrong) then added to 33 to make a number in the Ascii printable character range. As you can see from the table below, alphabetical letters - good, numbers – ok, most special characters – bad (except :;<=>?@).

See the Wikipedia FASTQ format page for more information.

What character in the quality score string in the fastq entry above represents the best base quality?

About compressed files

Sequencing data files can be very large - from a few megabytes to gigabytes. And with NGS giving us longer reads and deeper sequencing at decreasing price points, it's not hard to run out of storage space. As a result, most sequencing facilities will give you compressed sequencing data files. The most common compression program used for individual files is gzip and its counterpart gunzip whose compressed files have the .gz extension. The tar and zip programs are most commonly used for compressing directories.

Let's take a look at the size difference between uncompressed and compressed files. We use the -l option of ls to get a long listing that includes the file size, and -h to have that size displayed in "human readable" form rather than in raw byte sizes.

ls -lh $BI/web/yeast_stuff/*.fastq
ls -lh $BI/web/yeast_stuff/*.fastq.gz

About how big are the compressed files? The uncompressed files? About what is the compression factor?

You may be tempted to want to uncompress your sequencing files in order to manipulate them more directly – but resist that temptation. Nearly all modern bioinformatics tools are able to work on .gz files, and there are tools and techniques for working with the contents of compressed files without ever uncompressing them.

gzip and gunzip

With no options, gzip compresses the file you give it in-place. Once all the content has been compressed, the original uncompressed file is removed, leaving only the compressed version (the original file name plus a .gz extension). The gzunzip function works in a similar manner, except that its input is a compressed file with a .gz file and produces an uncompressed file without the .gz extension.

# make sure you're in your $SCRATCH/core_ngs/fastq_prep directory
cp $CLASSDIR/common/small.fq .

# check the size, then compress it in-place
ls -lh 
gzip small.fq

# check the compressed file size, then uncompress it
ls -lh 
gunzip small.fq.gz

head and tail, more or less

One of the challenges of dealing with large data files, whether compressed or not, is finding your way around the data – finding and looking at relevant pieces of it. Except for the smallest of files, you can't open them up in a text editor because those programs read the whole file into memory, so will choke on sequencing data files! Instead we use various techniques to look at pieces of the files at a time.

The first technique is the use of "pages" – we've already seen this with the more command. Review its use now on our small uncompressed file:

# Use spacebar to advance a page; Ctrl-c to exit
more small.fz 

Another pager, with additional features, is less. The most useful feature of less is the ability to search – but it still doesn't load the whole file into memory, so searching a really big file can be slow.

Here's a summary of the most common less navigation commands, once the less pager is active. It has tons of other options (try less --help).

• q – quit
• Ctrl-f or space – page forward
• Ctrl-b – page backward
• /<pattern> – search for <pattern> in forward direction
  • n – next match
  • N – previous match
• ?<pattern> – search for <pattern> in backward direction
  • n – previous match going back
  • N – next match going forward

If you start less with the -N option, it will display line numbers.

What line of small.fq contains the read name with grid coordinates 2316:10009:100563?

less -N small.fq
/2316:10009:100563

less -N small.fq
/2316:10009:100563

line number 905, which looks like this:
905 @HWI-ST1097:127:C0W5VACXX:5:2316:10009:100563 1:N:0:CTCAGA

For a really quick peak at the first few lines of your data, there's nothing like the head command. By default it displays the first 10 lines of data from the file you give it or from its standard input. With an argument -NNN (that is a dash followed by some number), it will show that many lines of data.

# shows 1st 10 lines
head small.fq

# shows 1st 100 lines -- might want to pipe this to more to see a bit at a time
head -100 small.fq | more

The yang to head's ying is tail, which by default it displays the last 10 lines of its data, and also uses the -NNN syntax to show the last NNN lines. (Note that with very large files it may take a while for tail to start producing output because it has to read through the file sequentially to get to the end.)

But what's really cool about tail is its -n +NNN syntax. This displays all the lines starting at line NNN. Note this syntax: the -n option switch follows by a plus sign ( + ) in front of a number – the plus sign is what says "starting at this line"! Try these examples:

# shows the last 10 lines
tail small.fq

# shows the last 100 lines -- might want to pipe this to more to see a bit at a time
tail -100 small.fq | more

# shows all the lines starting at line 900 -- better pipe it to more!
tail -n +900 small.fq | more

# shows 15 lines starting at line 900 because we pipe to head -15
tail -n +900 small.fq | head -15

gunzip -c tricks

Ok, now you know how to navigate an uncompressed file using head and tail, more or less. But what if your FASTQ file has been compressed by gzip? You don't want to uncompress the file, remember? So you use the gunzip -c trick. This uncompresses the file, but instead of writing the uncompressed data to another file (without the .gz extension) it write it to its standard output where it can be piped to programs like your friends head and tail, more or less.

Let's illustrate this using one of the compressed files in your fq subdirectory:

gunzip -c fq/Sample_Yeast_L005_R1.cat.fastq.gz | more
gunzip -c fq/Sample_Yeast_L005_R1.cat.fastq.gz | head
gunzip -c fq/Sample_Yeast_L005_R1.cat.fastq.gz | tail
gunzip -c fq/Sample_Yeast_L005_R1.cat.fastq.gz | tail -n +900 | head -15

Counting your sequences

One of the first thing to check is that your FASTQ files are the same length, and that length is evenly divisible by 4. The wc command (word count) using the -l switch to tell it to count lines, not words, is perfect for this. It's so handy that you'll end up using wc -l a lot to count things. It's especially powerful when used with filename wildcarding.

wc -l small.fq
head -100 small.fq > small2.fq
wc -l small*.fq

You can also pipe the output of gunzip -c to wc -l to count lines in your compressed FASTQ file.

How many lines are in the Sample_Yeast_L005_R1.cat.fastq.gz file? How many sequences is this?

gunzip -c fq/Sample_Yeast_L005_R1.cat.fastq.gz | wc -l

How to do math on the command line

The bash shell has a really strange syntax for arithmetic: it uses a double-parenthesis operator after the '$' sign (which means evaluate this expression). Go figure.

echo $((2368720 / 4))

Remember how the backticks evaluate the enclosed expression and substitute its output into a string? Here's how you would combine this math expression with gunzip -c line counting on your file using the magic of backtick evaluation:

echo "$((`gunzip -c fq/Sample_Yeast_L005_R1.cat.fastq.gz | wc -l` / 4))"

Whew!

Processing multiple compressed files

You've probably figured out by now that you can't easily use filename wildcarding along with gunzip -c and piping to process multiple files. For this, you need to code a for loop in bash. Fortunately, this is pretty easy. Try this:

for fname in fq/*.gz; do
   echo "$fname has $((`gunzip -c $fname | wc -l` / 4)) sequences"
done

Here fname is the name I gave the variable that is assigned a different file generated by the filename wildcarding list, each time through the loop. The actual file is then referenced as $file inside the loop.

Note the general structure of the for loop. Different portions of the structure can be separated on different lines (like <something> and <something else> below) or put on one line separated with a semicolon ( ; ) like before the do keyword below.

for <variable name> in <expression>; do 
  <something>
  <something else>
done

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

# setup
cds
mkdir samstat_test
cd samstat_test
cp $BI/ngs_course/intro_to_mapping/data/SRR030257_1.fastq .

# run the program
$BI/bin/samstat SRR030257_1.fastq

This produces a file named SRR030257_1.fastq.html which you need to view in a web browser. We put a copy at this URL:

http://loving.corral.tacc.utexas.edu/bioiteam/SRR030257_1.fastq.html

Running samstat on BAM files

# setup
cds
mkdir samstat_test2
cd samstat_test2
cp $BI/web/yeast_stuff/yeast_chip_sort.bam .

# run the program
$BI/bin/samstat yeast_chip_sort.bam

This produces a file named yeast_chip_sort.bam.html which you need to view in a web browser. We put a copy at this URL:

http://loving.corral.tacc.utexas.edu/bioiteam/yeast_stuff/yeast_chip_sort.bam.html

Note that by default, samstat only considers mapped reads for BAM files, although this behavior can be changed by piping the subset of reads you want analyzed to samstat from a samtools view command.

FASTQ Manipulation 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.

FASTX-Toolkit is available via the TACC module system.

module spider fastx_toolkit
module load fastx_toolkit

Here's an example of how to run fastx_trimmer to trim all input sequences down to 50 bases. By default the program reads its input data from standard input and writes trimmed sequences to standard output:

gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz | fastx_trimmer -l 50 -Q 33 > trimmed.fq

• The -l 50 option says that base 50 should be the last base (i.e., trim down to 50 bases)
• the -Q 33 option specifies how base qualities on the 4th line of each fastq entry are encoded. The FASTX toolkit is an older program, written in the time when Illumina base qualities were encoded differently. These days Illumina base qualities follow the Sanger FASTQ standard (Phred score + 33 to make an ASCII character).

Exercise: compressing the fastx_trimmer output

How would you tell fastx_trimmer to compress (gzip) its output file?

gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz | fastx_trimmer -l 50 -Q 33 -z > trimmed.fq.gz

gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz | fastx_trimmer -l 50 -Q 33 | gzip > trimmed.fq.gz

Exercise: fastx toolkit programs

What other fastx manipulation programs are part of the fastx toolkit?

ls $TACC_FASTX_BIN

Adapter trimming

Data from RNA-seq or other library prep methods that resulted in very short fragments can cause problems with moderately long (50-100bp) reads since the 3' end of sequence can be read through to the 3' adapter at a variable position. This 3' adapter contamination can cause the "reql" insert sequence not to align because the adapter sequence does not correspond to the bases at the 3' end of the reference genome sequence.

Unlike general fixed-length trimming (e.g. trimming 100 bp sequences to 40 or 50 bp), adapter trimming removes differing numbers of 3' bases depending on where the adapter sequence is found.

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.

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.

cutadapt -m 22 -O 10 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

cutadapt -m 22 -O 10 -a TGATCGTCGGACTGTAGAACTCTGAACGTGTAGA

Notes:

• The -m 22 option says to discard any sequence that is smaller than 22 bases after trimming. This avoids problems trying to map very short, highly ambiguous sequences.
• the -O 10 option says not to trim 3' adapter sequences unless at least the first 10 bases of the adapter are ssen at the 3' end of the read. This prevents trimming short 3' sequences that just happen by chance to match the first few adapter sequence bases.

<P5 capture> <indexRead2> <Read 1 primer> [insert] <Read 2 primer> <indexRead1> <P7 capture>

AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

TCTACACGTTCAGAGTTCTACAGTCCGACGATCA    # small RNA sequencing primer site
CAGGTTCAGAGTTCTACAGTCCGACGATCA        # "other"
TCTACACTCTTTCCCTACACGACGCTCTTCCGATCT  # TruSeq Read 1 primer site. This is the RC of the R2 adapter

raries, both R1 and R2 reads}
cutadapt -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC

TCTACACGTTCAGAGTTCTACAGTCCGACGATCA    # R1 primer - small RNA sequencing Read 1 primer site, 5' to 3'
TGATCGTCGGACTGTAGAACTCTGAACGTGTAGA    # R2 adapter contaminent (RC of R1 small RNA sequencing Read 1 primer)

Flexbar

Flexbar provides a flexible suite of commands for demultiplexing barcoded reads and removing adapter sequences from the ends of reads.

flexbar -n 1 --adapters adaptors.fna --source example.fastq --target example.ar --format fastq-sanger --adapter-threshold 2 --adapter-min-overlap 6 --adapter-trim-end RIGHT_TAIL

>adaptor1
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
>adaptor2
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTG
>adaptor1_RC
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
>adaptor2_RC
CAAGCAGAAGACGGCATACGAGATNNNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT

Note that flexbar only searches for the exact 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 with the potential benefit that many illumina adaptor sequences are already "built-in". It is available here.