FASTQ Quality Control Tutorial

Download the datasets

Datasets can be downloaded from: https://zenodo.org/records/21871083

Structure of FASTQ files

FASTQ files are the output from sequencing platforms (like Illumina) and typically consist of four lines per read:

@SRR15369215.126490887         # Sequence identifier
GGACCTTCTGTCA...              # Nucleotide sequence
+                             # Separator
AAFFFJJJJJJJJ...              # Base call quality scores

1. Count number of reads

For SRR10764405_1_subset.fastq.gz:

zcat sub_samples2/SRR10764405_1_subset.fastq.gz | wc -l | awk '{print $1 / 4}'

For SRR10764405_2_subset.fastq.gz:

zcat sub_samples2/SRR10764405_2_subset.fastq.gz | wc -l | awk '{print $1 / 4}'

For all FASTQ files:

# Strategy 1: count the total number of lines and divide by 4
# since all FASTQ reads are represented by 4 lines

for file in sub_samples2/*.fastq.gz; do
  read_count=$(( $(zcat "$file" | wc -l) / 4 ))
  echo "$(basename "$file"): $read_count reads"
done

# Strategy 2: count the number of lines with +
# since every FASTQ read has one line containing only +

zcat sub_samples2/SRR10764405_1_subset.fastq.gz | grep "^+$" | wc -l

2. Reads shorter than 100 bp

For SRR10764405_1_subset.fastq.gz:

zcat sub_samples2/SRR10764405_1_subset.fastq.gz | awk 'NR % 4 == 2 {print length($0)}' | awk '$1<100' | wc -l

For SRR10764405_2_subset.fastq.gz:

zcat sub_samples2/SRR10764405_2_subset.fastq.gz | awk 'NR % 4 == 2 {print length($0)}' | awk '$1<100' | wc -l

For all files:

for file in sub_samples2/*.fastq.gz; do
  reads=$(zcat "$file" | awk 'NR % 4 == 2 {print length($0)}' | awk '$1<100' | wc -l)
  echo "$(basename "$file"): $reads reads shorter than 100bp"
done

3. Reads longer than 100 bp

For SRR10764405_1_subset.fastq.gz:

zcat sub_samples2/SRR10764405_1_subset.fastq.gz | awk 'NR % 4 == 2 {print length($0)}' | awk '$1>100' | wc -l

For SRR10764405_2_subset.fastq.gz:

zcat sub_samples2/SRR10764405_2_subset.fastq.gz | awk 'NR % 4 == 2 {print length($0)}' | awk '$1>100' | wc -l

For all files:

for file in sub_samples2/*.fastq.gz; do
  reads=$(zcat "$file" | awk 'NR % 4 == 2 {print length($0)}' | awk '$1>100' | wc -l)
  echo "$(basename "$file"): $reads reads longer than 100bp"
done

4. Reads not equal to 100 bp

For SRR10764405_1_subset.fastq.gz:

zcat sub_samples2/SRR10764405_1_subset.fastq.gz | awk 'NR % 4 == 2 {if (length($0) != 100) print length($0)}' | wc -l

For SRR10764405_2_subset.fastq.gz:

zcat sub_samples2/SRR10764405_2_subset.fastq.gz | awk 'NR % 4 == 2 {if (length($0) != 100) print length($0)}' | wc -l

For all files:

for file in sub_samples2/*.fastq.gz; do
  count=$(zcat "$file" | awk 'NR % 4 == 2 {if (length($0) != 100) print length($0)}' | wc -l)
  echo "$(basename "$file"): $count reads not equal to 100bp"
done

5. Reads contaminated with Illumina adapters

Adapter sequence: CTGTCTCTTATACACATCT

For SRR10764405_1_subset.fastq.gz:

zcat sub_samples2/SRR10764405_1_subset.fastq.gz | awk 'NR % 4 == 2' | grep -c 'CTGTCTCTTATACACATCT'

For SRR10764405_2_subset.fastq.gz:

zcat sub_samples2/SRR10764405_2_subset.fastq.gz | awk 'NR % 4 == 2' | grep -c 'CTGTCTCTTATACACATCT'

Note that these commands look for a complete match to the adapter sequence. Most reads will contain only a partial length of the adapter sequence.

For all files:

for file in sub_samples2/*.fastq.gz; do
  count=$(zcat "$file" | awk 'NR % 4 == 2' | grep -c 'CTGTCTCTTATACACATCT')
  echo "$file: $count"
done

Expected: All files should report 0 contamination for complete adapter matches.