Page Comparison

...

What the 2nd sequence in the file $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_1.fastq is?

head $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_1.fastq 

...

wc -l $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_1.fastq 

Exercise: Counting FASTQ file lines

...

gunzip -c $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_2.fastq.gz | wc -l

...

mkdir $SCRATCH/BDIB_fastqc_tutorial
cd $SCRATCH/BDIB_fastqc_tutorial
cp $WORK$BI/GVA2016gva_course/mapping/data/SRR030257_1.fastq .
module load fastqc
 
fastqc -h  # examine program options
fastqc SRR030257_1.fastq  # run the program

...

...

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

...

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.

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.

cat SRR030257_1.fastq | fastx_trimmer -l 34 -Q 33 > 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 -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).
• 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

How would you tell fastx_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?

cat SRR030257_1.fastq | fastx_trimmer -l 34 -Q 

30 -z > SRR030257_1.trimmed.fastq.gz

cat SRR030257_1.fastq | fastx_trimmer -l 

34 -Q 

30 | gzip > SRR030257_1.trimmed.fastq.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 (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.

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)

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

$BI/gva_course/read_processing/JJM104_TAAGGCGA-TAGATCGC_L001_R1_001.fastq.gz
$BI/gva_course/read_processing/REL606.fna

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

$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

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?

fastx_clipper -i SRR030257_1.trimmed.fastq -o SRR030257_1.trimmed.depleted.fastq -a AAAAAAAAAAAAAAAAAAAA -l 34 -n 

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.

ls $TACC_FASTX_BIN

Consider how you might use fastq_quality_filter and fastq_quality_trimmer to limit your data based on quality scores.

Return to GVA2016 course page