2020 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 GSAF's Illumina HiSeq sequencer.

1 Setup
- 1.1 Set up directory for working with FASTQs
2 FASTQ Quality Assurance tools
- 2.1 FastQC
  - 2.1.1 Running FastQC
    - 2.1.1.1 Running fastqc on a FASTQ file
  - 2.1.2 Looking at FastQC output
- 2.2 Using MultiQC to consolidate multiple QC reports
- 2.3 FASTQ or BAM statistics with samstat
  - 2.3.1 Running samstat on FASTQ files
    - 2.3.1.1 Running samstat on FASTQ example
  - 2.3.2 Running samstat on BAM files
    - 2.3.2.1 Running samstat on BAM example
3 Trimming sequences
- 3.1 FASTX Toolkit
  - 3.1.1 Trimming FASTQ sequences to 50 bases with fastx_trimmer
- 3.2 Adapter trimming with cutadapt
  - - 3.2.1.1 cutadapt command for R1 sequences (GSAF RNA library)
    - 3.2.1.2 cutadapt command for R2 sequences (GSAF RNA library)
    - 3.2.1.3 Illumina library read layout
    - 3.2.1.4 Read 2 primer, 5' to 3', used as R1 sequence adapter
    - 3.2.1.5 Read 1 primer depends on library construction
    - 3.2.1.6 Cutadapt adapter sequence for ChIP-seq lib
    - 3.2.1.7 Small RNA library Read 1 primer, 5' to 3', used as R2 sequence adapter
  - 3.2.2 cutadapt example
    - 3.2.2.1 Setup for cutadapt on miRNA FASTQ
    - 3.2.2.2 Cutadapt batch command for R1 FASTQ
    - 3.2.2.3 cutadapt log file
  - 3.2.3 paired-end data considerations
  - 3.2.4 running cutadapt in a batch job
    - 3.2.4.1 Setup for cutadapt
    - 3.2.4.2 trim_adapters.sh
    - 3.2.4.3 Adapter trimming commands file
    - 3.2.4.4 Create and submit cutadapt batch script
    - 3.2.4.5 cutadapt log file
    - 3.2.4.6 cutadapt log file

Setup

Make sure you are logged in to our dedicated login node: login5.ls5.tacc.utexas.edu.

Set up to process the yeast data if you haven't already.

Set up directory for working with FASTQs

# Create a $SCRATCH area to work on data for this course,
# with a sub-direct[1ory for pre-processing raw fastq files
mkdir -p $SCRATCH/core_ngs/fastq_prep

# Make a symbolic links to the original yeast data:
cd $SCRATCH/core_ngs/fastq_prep
ln -s -f /work/projects/BioITeam/projects/courses/Core_NGS_Tools/yeast_stuff/Sample_Yeast_L005_R1.cat.fastq.gz
ln -s -f /work/projects/BioITeam/projects/courses/Core_NGS_Tools/yeast_stuff/Sample_Yeast_L005_R2.cat.fastq.gz

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:

Should I trim low quality bases?
- consult the Per base sequence quality report
  - based on all sequences
Do I need to remove adapter sequences?
- consult the Adapter Content report
Do I have other contamination?
- consult the Overrepresented Sequences report
  - based on the 1st 100,000 sequences, trimmed to 75bp
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

For many of its reports, FastQC analyzes only the first ~100,000 sequences in order to keep processing and memory requirements down. Consult the Online documentation for each FastQC report for full details.

Running FastQC

FastQC is available in the TACC module system on ls5. To make it available:

module load fastqc

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

# make sure you're in your $SCRATCH/core_ngs/fastq_prep directory
cds
core_ngs/fastq_prep
fastqc small.fq

Exercise: What did FastQC create?

ls -l shows two new items. -rw-rw-r-- 1 abattenh G-801021 319652 Jun 13 16:50 small_fastqc.html -rw-rw-r-- 1 abattenh G-801021 344099 Jun 13 16:50 small_fastqc.zip small_fastqc.html is the FastQC report, in HTML format. small_fastqc.zip is a zipped (compressed) directory of FastQC output files.

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

unzip small_fastqc.zip

What was created?

ls -l shows one new item, the small_fastqc directory (note the "d" in "drwxrwxr-x") drwxrwxr-x 4 abattenh G-801021 4096 Jun 13 16:50 small_fastqc ls -l small_fastqc shows the directory contents: drwxrwxr-x 2 abattenh G-801021 4096 Jun 13 16:50 Icons drwxrwxr-x 2 abattenh G-801021 4096 Jun 13 16:50 Images -rw-rw-r-- 1 abattenh G-801021 77561 Jun 13 16:50 fastqc.fo -rw-rw-r-- 1 abattenh G-801021 25635 Jun 13 16:50 fastqc_data.txt -rw-rw-r-- 1 abattenh G-801021 319652 Jun 13 16:50 fastqc_report.html -rw-rw-r-- 1 abattenh G-801021 446 Jun 13 16:50 summary.txt

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 laptap (read more at Copying files from TACC to your laptop).

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

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

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 paried-end ATAC-seq datasets – 4 FASTQ files total. First stage the data:

mkdir -p $SCRATCH/core_ngs/multiqc/fqc.atacseq
cd $SCRATCH/core_ngs/multiqc/fqc.atacseq
cp $CORENGS/multiqc/fqc.atacseq/*.zip .

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

50knuclei_S56_L007_R1_001_fastqc.zip  5knuclei_S77_L008_R1_001_fastqc.zip  
50knuclei_S56_L007_R2_001_fastqc.zip  5knuclei_S77_L008_R2_001_fastqc.zip

Now make the MultiQC accessible in your environment. It is not in the standard TACC module system, but is in the BioContainers.

# Load the main BioContainers module, then the multiqc module
module load biocontainers  # may take a while
module load multiqc

# Ask multiqc for its usage information
multiqc --help

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:

cd $SCRATCH/core_ngs/multiqc
multiqc fqc.atacseq

Exercise: How many reports did multiqc find?

Based on its execution output, it found 4 reports [WARNING] multiqc : MultiQC Version v1.9 now available! [INFO ] multiqc : This is MultiQC v1.7 [INFO ] multiqc : Template : default [INFO ] multiqc : Searching 'fqc.atacseq/' [INFO ] fastqc : Found 4 reports [INFO ] multiqc : Compressing plot data [INFO ] multiqc : Report : multiqc_report.html [INFO ] multiqc : Data : multiqc_data [INFO ] multiqc : MultiQC complete

Exercise: What was created by running multiqc?

One file was created (multiqc_report.html) and one directory (multiqc_data).