Pre-processing raw 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 one of GSAF's Illumina sequencers.

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, and why:

1. Should I trim low quality bases?

   • consult the Per base sequence quality report

     • based on all sequences

2. Do I need to remove adapter sequences?

   • consult the Adapter Content report

3. Do I have other contamination?
   
   

   • consult the Overrepresented Sequences report

     • based on the 1st 100,000 sequences, trimmed to 75bp

4. How complex is my library?

   • consult the Sequence Duplication Levels report

   • but remember that different experiment types are expected to have vastly different duplication profiles
     
     
     

Running FastQC

Make sure you're in an idev session. If you're in an idev session, the hostname command will display a name like c455-021.ls6.tacc.utexas.edu. But if you're on a login node the hostname will be something like login2.ls6.tacc.utexas.edu.

If you're on a login node, start an idev session like this:

Start an idev session

FastQC is available as part of BioContainers on ls6. To make it available:

It has a number of options (see fastqc --help | more) but can be run very simply with just a FASTQ file as its argument.

Running fastqc on a FASTQ file

Exercise: What did FastQC create?

Let's unzip the .zip file and see what's in it.

What was created?

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. One way to do this is to copy the results back to your laptop, for example by using scp from your computer (read more at Copying files from TACC to your laptop).

For convenience, we put an example FastQC report at this URL:
https://web.corral.tacc.utexas.edu/BioinformaticsResource/CoreNGS/yeast_stuff/Sample_Yeast_L005_R1.cat_fastqc/fastqc_report.html 

Exercise: Based on this FastQC output, should we trim this data?

Using MultiQC to consolidate multiple QC reports

FastQC reports are all well and good, but what if you have dozens of samples? It quickly becomes tedious to have to look through all the separate FastQC reports, including separate R1 and R2 reports for paired end datasets.

The MultiQC tool helps address this issue. Once FastQC reports have been generated, it can scan them and create a consolidated report from all the individual reports.

Whats even cooler, is that MultiQC can also consolidate reports from other bioinformatics tools (e.g. bowtie2 aligner statistics, samtools statistics, cutadapt, Picard, and may more). And if your favorite tool is not known by MultiQC, you can configure custom reports fairly easily. For more information, see this recent Byte Club tutorial on Using MultiQC.

Here we're just going to create a MultiQC report for two paired-end ATAC-seq datasets – 4 FASTQ files total. First stage the data:

You should see these 4 files in your $SCRATCH/core_ngs/multiqc/fqc.atacseq directory:

Now make the BioContainers MultiQC accessible in your environment.

Even though multiqc has many options, it is quite easy to create a basic report by just pointing it to the directory where individual reports are located:

Exercise: How many reports did multiqc find?

Exercise: What was created by running multiqc?

You can see the resulting MultiQC report here: https://web.corral.tacc.utexas.edu/BioinformaticsResource/CoreNGS/reports/atacseq/multiqc_report.html.

And an example of a MultiQC report that includes both standard and custom plots is this is the Tag-Seq post-processing MultiQC report produced by the Bioinformatics Consulting Group: https://web.corral.tacc.utexas.edu/BioinformaticsResource/CoreNGS/reports/mqc_tagseq_trim_JA21030_SA21045_mouse.html

Trimming sequences

There are two main reasons you may want to trim your sequences:

• As a quick way to remove 3' adapter contamination, when extra bases provide little additional information

  • For example, 75+ bp ChIP-seq reads – 50 bases are more than enough for a good mapping, and trimming to 50 is easier than adapter removal, especially for paired end data.

  • You would not choose this approach for RNA-seq data, where 3' bases may map to a different exon, and that is valuable information.

    • Instead you would specifically remove adapter sequences.

• Low quality base reads from the sequencer can affect some programs
  
  

  • This is an issue with sequencing for genome or transcriptome assembly.

  • Aligners such as bwa and bowtie2 seem to do fine with a few low quality bases, soft clipping them if necessary.

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 both 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 as a BioContainers module.

Here's an example of how to run fastx_trimmer to trim all input sequences down to 50 bases.

Where does fastx_trimmer read its input from? And where does it write its output? Ask the program for its usage.

The usage: is its help information

Because the [-i INFILE] [-o OUTFILE] options are shown in brackets [ ], reading from a file and writing to a file are optional. That means that by default the program reads its input data from standard input and writes trimmed sequences to standard output:

Trimming FASTQ sequences to 50 bases with fastx_trimmer

• 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, so its default does not work for modern FASTQ files.

    • These days Illumina base qualities follow the Sanger FASTQ standard (Phred score + 33 to make an ASCII character).

    • This option is not really required here because we're just hard trimming, so the program doesn't have to interpret the quality scores. But the -Q 33 option would be required if you were trimming according to base qualities.

  • Note that the fastq_trimmer help does not document this -Q option!

Exercise: compressing fastx_trimmer output

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

Exercise: other fastx toolkit programs

What other FASTQ manipulation programs are part of the FASTX Toolkit?

The FASTX Toolkit also has programs that work on FASTA files. To see them, type fasta_ then tab twice (completion) to see their names.

Adapter trimming with cutadapt

Data from RNA-seq or other library prep methods that result in short fragments can cause problems with moderately long (50-100bp) reads, since the 3' end of sequences can be read into (or even through) to the 3' adapter at different read offsets . 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 50 bp), specific adapter trimming removes differing numbers of 3' bases depending on where the adapter sequence is found.

The cutadapt program, available in BioContainers, is an excellent tool for removing adapter contamination.

A common application of cutadapt is to remove adapter contamination from RNA library sequence data. Here we'll show that for some small RNA libraries sequenced by GSAF, using their documented small RNA library adapters.

When you run cutadapt you give it the adapter sequence to trim, and the adapter sequence is different for R1 and R2 reads. Here's what the options look like (without running it on our files yet).

cutadapt command for R1 sequences (GSAF RNA library)

cutadapt command for R2 sequences (GSAF RNA library)

Notes:

• The -m 22 option says to discard any sequence that is smaller than 22 bases (minimum) after trimming.
  
  

  • This avoids problems trying to map very short, highly ambiguous sequences.

• the -O 4 (Overlap) option says not to trim 3' adapter sequences unless at least the first 4 bases of the adapter are seen 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.

Figuring out which adapter sequence to use when can be tricky. Your sequencing provider can tell you what adapters they used to prep your libraries. For GSAF's adapter layout, please refer to Illumina - all flavors (USE with Caution, this is outdated but can be useful for a basic understanding of the adapters, the GSAF primarily only uses UDI's for all projects) (you may want to read all the "gory details" below later).

Exercise: other cutadapt options

The cutadapt program has many options. Let's explore a few.

How would you tell cutadapt to trim trailing N's?

How would you control the accuracy (error rate) of cutadapt's matching between the adapter sequences and the FASTQ sequences?

Suppose you are processing 100 bp reads with 30 bp adapters. By default, how many mismatches between the adapter and a sequence will be tolerated?

How would you require a more stringent matching (i.e., allowing fewer mismatches)?

cutadapt example

Let's run cutadapt on some real human miRNA (micro-RNA) data.

First, stage the data we want to use. This data is from a small RNA library where the expected insert size is around 15-25 bp.

Setup for cutadapt on miRNA FASTQ

Exercise: How many reads are in these files? Is it single end or paired end data?

Exercise: How long are the reads?

Adapter trimming is a rather slow process, and these are large files. So to start with we're going to create a smaller FASTQ file to work with.

Now execute cutadapt like this. Note that the backslash ( \ ) here is just a line continuation character so that we can split a long command onto multiple lines to make it more readable.

Cutadapt command for R1 FASTQ