...
Warning | |||||||
---|---|---|---|---|---|---|---|
As we discussed yesterday, we do not want to run things on the head node, so we should start a new idev session. The reservation name has changed from CCBB_Bio_Summer_School_2016_Day1 to CCBB_Bio_Summer_School_2016_Day2. See if you can figure out how to start your own idev session.
|
...
See the Wikipedia FASTQ format page for more information.
Now that you know the basics, see if you can complete the following exercises on your own.
Exercise: Examine the 2nd sequence in a FASTQ file
What is the 2nd sequence in the file $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_1.fastq
is?
Expand | |||||||
---|---|---|---|---|---|---|---|
| |||||||
Use the
|
Expand | ||
---|---|---|
| ||
The 2nd sequence has ID = If thats what you thought it was congratulations, if it is different, do you see where we got it from? If it doesn't make sense ask for help. |
...
Code Block | ||||
---|---|---|---|---|
| ||||
wc -l $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_1.fastq |
Exercise: Counting FASTQ file lines
...
Code Block | ||||
---|---|---|---|---|
| ||||
gunzip -c $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_2.fastq.gz | wc -l |
...
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.
FASTQ Evaluation Tools
The Place your sticky note on your computer when you have made it this far and start looking over the fastqc links below. Once everyone has caught up we will go over this together.
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 is a tool that produces a quality analysis report on FASTQ files.
Useful links:
...
Online documentation for FastQC
First and foremost, the FastQC "Summary" on the left 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.
...
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.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:. We don't want to clutter up our work space so copy the SRR030257_1.fastq file to a new directory named BDIB_fastqc_tutorial on scratch, use the module system to load fastqc, use fastqc's help option after the module is loaded 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)
Code Block | ||||||
---|---|---|---|---|---|---|
| # setup
module load fastqc
cds
mkdir fastqc_test
cd fastqc_test
| |||||
mkdir $SCRATCH/BDIB_fastqc_tutorial cd $SCRATCH/BDIB_fastqc_tutorial cp $BI/gva_course/webmapping/yeast_stuffdata/Sample_Yeast_L005_R1.catSRR030257_1.fastq .gz module .load fastqc # runningfastqc the-h program fastqc Sample_Yeast_L005_R1.cat.fastq.gz # examine program options fastqc SRR030257_1.fastq # examinerun extra options fastqc -hthe program |
Exercise: FastQC results
What did FastQC create?
Expand | |||||
---|---|---|---|---|---|
| |||||
The Sample_Yeast_L005_R1.catSRR030257_1.fastq .gz file is what we analyzed, so FastQC created the other two items. SampleSRR030257_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. |
...
1_fastqc.html represents the results in a file viewable in a web browser. SRR030257_1_fastqc.zip is just a Zipped (compressed) version of the results. |
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 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.If you want to skip that step (we recommend doing it for practice!), we have put a copy of the output at this URL:
Code Block | |||||||
---|---|---|---|---|---|---|---|
| http://web.corral
| ||||||
# on tacc terminal pwd # on new terminal of local computer scp <username>@ls5.tacc.utexas.edu/BioITeam/yeast_stuff/Sample_Yeast_L005_R1.cat_fastqc/fastqc_report.html :<pwd_results_from_other_window>/SRR030257_1_fastqc.html ~/Desktop # open the newly transfered file from from the desktop and see how the data looks |
Exercise: Should we trim this data?
...
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
|
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:
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.
FASTX-Toolkit is available via the TACC module system.
Code Block | ||
---|---|---|
| ||
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:
No Format | ||
---|---|---|
| ||
gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz great, but more importantly, nearly 1.5% of all the sequences are all A's according to the Overrepresented sequences. This is something that often comes up in miseq data that has shorter insert sizes than the overall read length. Next we'll start looking at how to trim our data before continuing. |
FASTQ Processing Tools
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 and fastx_clipper for trimming specific sequences (such as adapters).
FASTX-Toolkit is available via the TACC module system.
Code Block | ||
---|---|---|
| ||
module spider fastx_toolkit
module load fastx_toolkit |
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, but fastx_trimmer has the advantage of being a module on TACC and therefore the easiest to use. By default the program reads its input data from standard input and writes trimmed sequences to standard output, use what you know about piping and printing text to the screen to trip the fastq file to 34 bases.
Expand | ||
---|---|---|
| ||
Type fastx_trimmer -h to see program documentation. Look below the possible solution for more detailed information of what to focus on |
Code Block | ||||
---|---|---|---|---|
| ||||
cat SRR030257_1.fastq | fastx_trimmer -l 34 > SRR030257_1.trimmed.fastq |
- The cat command prints the contents of the SRR030257_1.fastq file to the screen
- The | redirects that output to the fastx_trimmer command
- The -l 34 option says that base 34 should be the last base (i.e., trim down to 34 bases)
- The > redirects that output to the new file on the right side, in this case SRR030257_1.trimmed.fastq
Exercise: compressing the fastx_trimmer output
Compressed files are smaller, easier to transfer, and many programs allow for their use directly. How would you tell fastx_trimmer to compress (gzip) its output file?
Expand | ||
---|---|---|
| ||
Type fastx_trimmer -h to see program documentation |
Code Block | ||||
---|---|---|---|---|
| ||||
cat SRR030257_1.fastq | fastx_trimmer -l 34 -z > SRR030257_1.trimmed.fastq.gz |
Code Block | ||||
---|---|---|---|---|
| ||||
cat SRR030257_1.fastq | fastx_trimmer -l 5034 -Q| 33gzip > SRR030257_1.trimmed.fastq.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
...
gz |
Both of the above solutions give the same final product, but are clearly achieved in different ways. This is done to show you that 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 reproducable.
Adapter trimming
As mentioned above, fastx_clipper can be used to trim specific sequences, and based on our fastqc analysis, the sequence AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA is significantly overrepresented in our data. How would you use fastx_clipper to remove those sequences from the fastq file?
Expand | ||
---|---|---|
| ||
Type fastx_ trimmerclipper -h to see program documentation |
Expand | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|
| ||||||||||
You could supply the -z option like this:
Or you could gzip the output yourself:
|
Exercise: fastx toolkit programs
...
Look below the possible solution for more detailed information on what to focus on. |
Code Block | ||||
---|---|---|---|---|
| ||||
fastx_clipper -i SRR030257_1.trimmed.fastq -o SRR030257_1.trimmed.depleted.fastq -a AAAAAAAAAAAAAAAAAAAA -l 34 -n |
Command portion | purpose |
---|---|
-i SRR030257_1.trimmed.fastq | use this file as input |
-o SRR030257_1.trimmed.depleted.fastq | create this new output file |
-a AAAAAAAAAAAAAAAAA | remove bases containing this sequence |
-l 34 | discard any read shorter than 34 bases after sequence removed |
-n | keep reads containing "N" bases in them. Consider how they are treated in downstream applications |
Other fastx toolkit programs
What other fastx and fastq manipulation programs are part of the fastx toolkit?
...
The available modules are described on their website, but we can also learn some things from the command prompt using module spider fastx_toolkit.
Expand | |||||
---|---|---|---|---|---|
| |||||
Type fastx_ fast then tab twice 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.
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. 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.
Code Block | ||
---|---|---|
| ||
cutadapt -m 22 -O 10 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
|
Code Block | ||
---|---|---|
| ||
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.
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.
No Format | ||
---|---|---|
| ||
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
|
Code Block | ||
---|---|---|
| ||
>adaptor1
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
>adaptor2
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTG
>adaptor1_RC
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
>adaptor2_RC
CAAGCAGAAGACGGCATACGAGATNNNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
|
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.
More Example Data
See if you can figure out what's wrong with these data sets (copy them to your $SCRATCH
directory before analyzing them) and then process them to get rid of the problem(s). If you're very ambitious, you could also map them to the reference genomes and perform variant calling before and after cleaning them up to see how the results change. Each file has a different problem.
Example #1: Single-end Illumina MiSeq data for E. coli
Code Block | ||||
---|---|---|---|---|
| ||||
$BI/gva_course/read_processing/JJM104_TAAGGCGA-TAGATCGC_L001_R1_001.fastq.gz
$BI/gva_course/read_processing/REL606.fna |
Expand | ||
---|---|---|
| ||
This |
Example #2: Paired-end Illumina Genome Analyzer IIx data for E. coli
Code Block | ||||
---|---|---|---|---|
| ||||
$BI/gva_course/read_processing/61FTVAAXX_2_R1_ZDB172.fastq.gz
$BI/gva_course/read_processing/61FTVAAXX_2_R2_ZDB172.fastq.gz
$BI/gva_course/read_processing/REL606.fna |
...
title | What's wrong with this data? |
---|
...
Consider how you might use fastq_quality_filter and fastq_quality_trimmer to limit your data based on quality scores.
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.