Versions Compared

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

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/GVA2016/mapping/data/SRR030257_1.fastq is?

Expand
titleStuck? click here for a hint

Use the head command.

Code Block
titleStill stuck? Click here for the full Head command
collapsetrue
head $WORK/GVA2016/mapping/data/SRR030257_1.fastq

...

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. 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:

...

titleRunning FastQC example

...

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
titleRunning FastQC example
collapsetrue
mkdir $SCRATCH/BDIB_fastqc_tutorial
cd $SCRATCH/BDIB_fastqc_tutorial
cp $WORK/GVA2016/mapping/data/SRR030257_1.fastq .
module load fastqc
 
fastqc -h  # examine program options
fastqc SRR030257_1.fastq  # run the program

Exercise: FastQC results

What did FastQC create?

Expand
titleAnswer
No Format
titlels -l shows something like this
drwxrwxr-rwxr-xr-x 41 abattenhded G-803889802740 498588268    4096 May 2023 2212:5906 Sample_Yeast_L005_R1.cat_fastqcSRR030257_1.fastq
-rw-r-rw-r-- 1 abattenhded G-803889802740    198239291714 May 2023 2212:5907 SampleSRR030257_Yeast_L005_R1.cat_1_fastqc.ziphtml
-rwxr-xr-x-rw-r--r-- 1 abattenhded G-803889802740 51065629 May 20 22:59 Sample_Yeast_L005_R1.cat.fastq.gz
The Sample_Yeast_L005_R1.cat.fastq.gz
 455677 May 23 12:07 SRR030257_1_fastqc.zip

The SRR030257_1.fastq 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 directory1_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 the results back to your 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:

http://web.corral
Code Block
titleFastQC results URL
Transferring fastqc data back to computer
collapsetrue
# 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
titleAnswer

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
titleRunning samstat on FASTQ example
# 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
titleURL for viewing samstat results
http://loving.corral.tacc.utexas.edu/bioiteam/SRR030257_1.fastq.html

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.

Additional FASTQ Processing Tools

Trimming low quality bases

...

No Format
titlefastx_trimmer example
gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gzcat SRR030257_1.fastq | fastx_trimmer -l 5034 -Q 33 > SRR030257_1.trimmed.fqfastq
  • The -l 5034 option says that base 50 34 should be the last base (i.e., trim down to 50 34 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).

...

Expand
titleAnswer

You could supply the -z option like this:

No Format
titlefastx_trimmer example
gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gzcat SRR030257_1.fastq | fastx_trimmer -l 5034 -Q 33 -z > SRR030257_1.trimmed.fqfastq.gz

Or you could gzip the output yourself:

gunzip -c $BI/web/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz
No Format
titlefastx_trimmer example
titlefastx_trimmer example
cat SRR030257_1.fastq | fastx_trimmer -l 50 -Q 33 | gzip > SRR030257_1.trimmed.fqfastq.gz

Exercise: fastx toolkit programs

...

Expand
titleThe gory details on the *-a* adapter sequence argument

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.

No Format
titleIllumina library read layout
<P5 capture> <indexRead2> <Read 1 primer> [insert] <Read 2 primer> <indexRead1> <P7 capture>

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:

Code Block
titleRead 2 primer, 5' to 3', used as R1 sequence adapter
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

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:

No Format
titleRead 1 primer depends on library construction
TCTACACGTTCAGAGTTCTACAGTCCGACGATCA    # small RNA sequencing primer site
CAGGTTCAGAGTTCTACAGTCCGACGATCA        # "other"
TCTACACTCTTTCCCTACACGACGCTCTTCCGATCT  # TruSeq Read 1 primer site. This is the RC of the R2 adapter

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.
Since it is the RC of the Read 2 primer, its RC is just the Read 1 primer back
Therefore, for ChIP-seq libraries only one cutadapt command is needed:

Code Block
titleCutadapt adapter sequence for ChIP-seq lib
raries, both R1 and R2 reads}
cutadapt -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC

For RNAseq libraries, we use the small RNA sequencing primer as the Read 1 primer.
The contaminent is then the RC of this, minus the 1st and last bases:

No Format
titleSmall RNA library Read 1 primer, 5' to 3', used as R2 sequence adapter
TCTACACGTTCAGAGTTCTACAGTCCGACGATCA    # R1 primer - small RNA sequencing Read 1 primer site, 5' to 3'
TGATCGTCGGACTGTAGAACTCTGAACGTGTAGA    # R2 adapter contaminent (RC of# R1 primer - small RNA sequencing 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
titleExample of trimming adaptor sequence from right end 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
Code Block
titleExample adaptors.fna file
>adaptor1
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
>adaptor2
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTG
>adaptor1_RC
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
>adaptor2_RC
CAAGCAGAAGACGGCATACGAGATNNNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT

...

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

Trimmomatic

Trimmomatic offers similar options to Flexbar with the potential benefit that many illumina adaptor sequences are already "built-in". It is available here.

...

Code Block
languagebash
titleExample read and reference files #1
$BI/gva_course/read_processing/JJM104_TAAGGCGA-TAGATCGC_L001_R1_001.fastq.gz
$BI/gva_course/read_processing/REL606.fna
Expand
titleWhat's wrong with this data?
This

 

Example #2: Paired-end Illumina Genome Analyzer IIx data for E. coli

...