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.
When following along here, please start an idev session for running any example commands:
|
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 |
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?
Use the head command. |
Executing the command above reports that the 2nd sequence has ID = @SRR030257.2 HWI-EAS_4_PE-FC20GCB:6:1:407:767/1, and the sequence TAAGCCAGTCGCCATGGAATATCTGCTTTATTTAGC |
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?
The wc -l command says there are 15200720 lines. FASTQ files have 4 lines per sequence, so the file has 15,200,720/4 or 3,800,180 sequences. |
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?
The wc -l command says there are 2368720 lines so the file has 2,368,720/4 or 592,180 sequences. |
The bash shell has a really strange syntax for arithmetic: it uses a double-parenthesis operator. Go figure.
|
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 is a tool that produces a quality analysis report on FASTQ files.
Useful links:
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:
|
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?
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. |
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?
The Per base sequence quality report does not look good. The data should probably be trimmed (to 40 or 50 bp) before alignment. |
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.
# 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 |
# 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.
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.
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 |
Exercise: compressing the fastx_trimmer output
How would you tell fastx_trimmer to compress (gzip) its output file?
Type fastx_trimmer -h to see program documentation |
You could supply the -z option like this:
Or you could gzip the output yourself:
|
Exercise: fastx toolkit programs
What other fastx manipulation programs are part of the fastx toolkit?
Type fastx_ then tab to see their names
|
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
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.
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:
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 "contaminent" 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 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.