Quality control of raw sequencing data, data cleaning, FASTQ file handling [GenomEast Wiki]

Quality control of raw sequencing data, data cleaning, FASTQ file handling

Instructor	Stephanie Le Gras
Duration	3.5 hours
Content	Quality control of raw sequencing data, data cleaning, FASTQ file handling (lecture)
Content	Practical session on quality controls of raw sequencing data (Hands-on)
Prerequisites	None

1 Prepare the analysis environment to use to perform the analysis

1.1 Create a working directory

cd /work/c-shd/[your login]
mkdir coursQC_Mapping

1.2 Move to the directory you've just created

cd coursQC_Mapping

1.3 Copy the datasets into your working directory

cp /user2/c-shd/sll13384/DU/Data/Sample_CRN-107_11/CRN-107_11-R*.fastq.gz .

1.4 Uncompress the files

gunzip CRN-107_11-R1.fastq.gz
gunzip CRN-107_11-R2.fastq.gz

1.5 Create a virtual working environment

mkdir Software
cd Software
wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh
bash Miniconda2-latest-Linux-x86_64.sh
# Press enter
# Press space until the end of the license
# Enter Yes
# Enter /work/c-shd/[your login]/Software/miniconda2
# Enter yes
conda config --add channels conda-forge
conda config --add channels defaults
conda config --add channels r
conda config --add channels bioconda
conda create --name dnaseq bedtools samtools bwa fastqc cutadapt picard gatk
# enter y
cp /user2/c-shd/sll13384/DU/Software/GenomeAnalysisTK-3.6.tar.bz2 .
source activate dnaseq
gatk-register $(pwd)/GenomeAnalysisTK-3.6.tar.bz2

1.6 Source working environment

source activate dnaseq
source /user2/c-shd/sll13384/DU/Scripts/envTD.sh

2 Compute the number of sequenced reads

File:

CRN-107_11-R1.fastq
CRN-107_11-R2.fastq

Tool:

Use the command line 'wc' ('wc' for word count)

Answer

wc -l CRN-107_11-R1.fastq CRN-107_11-R2.fastq

1122032 CRN-107_11-R1.fastq
1122032 CRN-107_11-R2.fastq

There is 1122032/4 = 280508 reads per file.

3 Assess the quality of the raw data

Files:

CRN-107_11-R1.fastq
CRN-107_11-R2.fastq

Tool:

fastqc

Tips

Put the results onto a directory called Fastqc

3.1 Create an output directory named Fastqc

Answer

mkdir Fastqc

3.2 Run the command fastqc on the two fastq files

By default, FastQC groups bases for reads >50bp. Disable this functionality.

Answer

fastqc --nogroup CRN-107_11-R1.fastq --outdir Fastqc
fastqc --nogroup CRN-107_11-R2.fastq --outdir Fastqc

3.3 Use firefox to look at the resulting html file (fastqc_report.html)

Answer

firefox Fastqc/CRN-107_11-R1_fastqc/fastqc_report.html
firefox Fastqc/CRN-107_11-R2_fastqc/fastqc_report.html

4 Trim the last nucleotide from sequences

Input files:

CRN-107_11-R1.fastq
CRN-107_11-R2.fastq

Output files:

CRN-107_11-R1_shorter.fastq
CRN-107_11-R2_shorter.fastq

Tool:

fastx_trimmer

Tips:

Be careful the encoding of the quality is Phred33 so tell it to fastx toolkit using the tag -Q 33.

Answer

fastx_trimmer -f 1 -l 100 -Q 33 -i CRN-107_11-R1.fastq -o CRN-107_11-R1_shorter.fastq
fastx_trimmer -f 1 -l 100 -Q 33 -i CRN-107_11-R2.fastq -o CRN-107_11-R2_shorter.fastq

4.1 Check the sequences length after trimming

Use a combination of 'head', 'tail' and 'wc'
Or open the file in a text editor with the column number displayed

Answer

head -2 CRN-107_11-R1_shorter.fastq | tail -1 | wc –c

head -2 CRN-107_11-R1.fastq | tail -1 | wc –c

5 Cut Illumina's adapters using Cutadapt

Input files:

CRN-107_11-R1_shorter.fastq
CRN-107_11-R2_shorter.fastq

Output files:

CRN-107_11-R1_trimmed.fastq
CRN-107_11-R2_trimmed.fastq

Tips:

The adaptater sequence is: AGATCGGAAGAGC
Sequence length after trimming should be higher that 30
Check out Cutadapt website for help on trimming adapters in paired end sequences.

5.1 Compute reverse complement of the adapter sequence using fastx_reverse_complement

Answer

echo -e ">illumina_adapter_forward\nAGATCGGAAGAGC" > adapterSeq.fa
fastx_reverse_complement -i adapterSeq.fa -o adapterSeqRevComp.fa

5.2 Run cutadapt

Answer

cutadapt -a AGATCGGAAGAGC -A GCTCTTCCGATCT -o CRN-107_11-R1_trimmed.fastq -p CRN-107_11-R2_trimmed.fastq CRN-107_11-R1_shorter.fastq CRN-107_11-R2_shorter.fastq

This is cutadapt 1.12 with Python 2.7.12
Command line parameters: -a AGATCGGAAGAGC -A GCTCTTCCGATCT -o CRN-107_11-R1_trimmed.fastq -p CRN-107_11-R2_trimmed.fastq CRN-107_11-R1_shorter.fastq CRN-107_11-R2_shorter.fastq
Trimming 2 adapters with at most 10.0% errors in paired-end mode ...
Finished in 18.71 s (67 us/read; 0.90 M reads/minute).

=== Summary ===

Total read pairs processed:            280,508
  Read 1 with adapter:                   7,798 (2.8%)
  Read 2 with adapter:                   4,267 (1.5%)
Pairs written (passing filters):       280,508 (100.0%)

Total basepairs processed:    56,101,600 bp
  Read 1:    28,050,800 bp
  Read 2:    28,050,800 bp
Total written (filtered):     56,061,352 bp (99.9%)
  Read 1:    28,024,747 bp
  Read 2:    28,036,605 bp

=== First read: Adapter 1 ===

Sequence: AGATCGGAAGAGC; Type: regular 3'; Length: 13; Trimmed: 7798 times.

No. of allowed errors:
0-9 bp: 0; 10-13 bp: 1

Bases preceding removed adapters:
  A: 31.7%
  C: 24.2%
  G: 23.0%
  T: 21.1%
  none/other: 0.0%

Overview of removed sequences
length	count	expect	max.err	error counts
3	6075	4382.9	0	6075
4	1425	1095.7	0	1425
5	252	273.9	0	252
6	20	68.5	0	20
8	2	4.3	0	2
9	2	1.1	0	0 2
10	9	0.3	1	1 8
11	2	0.1	1	2
12	2	0.0	1	0 2
35	1	0.0	1	1
38	1	0.0	1	0 1
40	1	0.0	1	0 1
53	1	0.0	1	0 1
58	1	0.0	1	0 1
66	1	0.0	1	0 1
88	1	0.0	1	0 1
100	2	0.0	1	0 2

=== Second read: Adapter 2 ===

Sequence: GCTCTTCCGATCT; Type: regular 3'; Length: 13; Trimmed: 4267 times.

No. of allowed errors:
0-9 bp: 0; 10-13 bp: 1

Bases preceding removed adapters:
  A: 31.4%
  C: 3.6%
  G: 25.6%
  T: 39.4%
  none/other: 0.0%

Overview of removed sequences
length	count	expect	max.err	error counts
3	3381	4382.9	0	3381
4	565	1095.7	0	565
5	250	273.9	0	250
6	54	68.5	0	54
7	12	17.1	0	12
9	2	1.1	0	0 2
10	1	0.3	1	0 1
44	1	0.0	1	0 1
62	1	0.0	1	0 1

6 Remove bad quality sequences

Files:

CRN-107_11-R1_trimmed.fastq
CRN-107_11-R2_trimmed.fastq

Tool:

DynamicTrim.pl (SolexaQA)

Tips:

Threshold for trimming: Phred score >= 10

6.1 Create an output directory for the data named SolexaQA

Answer

mkdir SolexaQA

6.2 Run the tool DynamicTrim.pl

Answer

DynamicTrim.pl -h 10 -d SolexaQA CRN-107_11-R1_trimmed.fastq CRN-107_11-R2_trimmed.fastq

Rename resulting files

mv SolexaQA/CRN-107_11-R1_trimmed.fastq.trimmed CRN-107_R1.fastq
mv SolexaQA/CRN-107_11-R2_trimmed.fastq.trimmed CRN-107_R2.fastq

7 Compress all output files

gzip CRN-107_11-R1.fastq CRN-107_11-R2.fastq
gzip CRN-107_11-R1_shorter.fastq CRN-107_11-R2_shorter.fastq
gzip CRN-107_11-R1_trimmed.fastq CRN-107_11-R2_trimmed.fastq
gzip CRN-107_R1.fastq CRN-107_R2.fastq

8 Put all temporary files into the directory named intermedFastqFiles

# Make output directory
mkdir intermedFastqFiles

# Move compressed files 
mv CRN-107_11-R1_shorter.fastq.gz CRN-107_11-R2_shorter.fastq.gz intermedFastqFiles 
mv CRN-107_11-R1_trimmed.fastq.gz CRN-107_11-R2_trimmed.fastq.gz intermedFastqFiles

9 Run FastQC on the final file

9.1 Create an output directory Fastqc_final

mkdir Fastqc_final

9.2 Run FastQC on the final fastq files

fastqc --nogroup CRN-107_R1.fastq.gz --outdir Fastqc_final