Quality Control with FastQC

Before you align or quantify anything, you inspect the raw reads. Sequencing introduces base-call errors, adapter contamination, and low-quality cycles, and feeding those artifacts downstream biases every result that follows. Quality control is the first checkpoint in any RNA-seq pipeline, and FastQC is the standard tool for it.

The FASTQ format

A FASTQ file stores each read as a record of exactly four lines. The first line starts with @ and holds the read identifier, the second line is the nucleotide sequence, the third line starts with + and optionally repeats the identifier, and the fourth line is a string of quality characters the same length as the sequence. Files are usually gzip-compressed, so peek at a record without decompressing the whole thing.

zcat raw_data/sample_R1.fastq.gz | head -n 4

@SRR1039508.1 HWI-ST177:1101:1225:2160/1
AGCTTGGCTACGGATCCGATTCGATCGTACGGATCA
+
IIIIIIIIIIIIIIIHHHHHHHGGGGGFFFFF===5

Line four encodes each base call's confidence as a Phred quality score, defined as Q = -10 * log10(P), where P is the probability that the base is wrong. Q30 means a 1 in 1000 error rate (99.9 percent accuracy) and Q20 is 1 in 100. Scores are stored as single ASCII characters using the Phred+33 offset, so a Q30 base prints as the character with ASCII code 63, a question mark, and the trailing 5 in the record above is a much lower Q20 base at the 3' end.

Running FastQC

mkdir -p qc_reports
fastqc -t 4 -o qc_reports raw_data/*.fastq.gz

Started analysis of sample_R1.fastq.gz
Started analysis of sample_R2.fastq.gz
Approx 25% complete for sample_R1.fastq.gz
Approx 25% complete for sample_R2.fastq.gz
Approx 50% complete for sample_R1.fastq.gz
Approx 50% complete for sample_R2.fastq.gz
Approx 75% complete for sample_R1.fastq.gz
Approx 75% complete for sample_R2.fastq.gz
Analysis complete for sample_R1.fastq.gz
Analysis complete for sample_R2.fastq.gz

For each input FastQC writes an HTML report plus a zip archive. The key modules are Per base sequence quality, whose boxplots should sit in the green band and typically dip at the 3' end; Per sequence GC content, which should be a single smooth peak, since a bimodal or shifted curve suggests contamination or rRNA; Adapter Content, which flags read-through adapter when fragments are shorter than the read length; and Overrepresented sequences, where you often see adapters, rRNA, or highly expressed transcripts. Treat the pass, warn, and fail icons as generic thresholds rather than verdicts: for RNA-seq a warn or even fail on Per base sequence content in the first ~12 bases and on Sequence Duplication Levels is expected, driven by random hexamer priming bias and highly expressed genes, not by bad data.

multiqc qc_reports/ -o multiqc_report