Quality Control of NGS Reads with FastQC and MultiQC
Chapter 4 of 5: MultiQC, automation and contamination
Step 5: Check again and combine everything with MultiQC
Run FastQC on the trimmed reads, then let MultiQC collect all reports into one page. MultiQC searches the folders you give it, recognises the output files of more than a hundred tools and puts their results side by side.
mkdir -p qc_trim
fastqc trim_R1.fastq.gz trim_R2.fastq.gz -o qc_trim/ -t 4
multiqc qc_raw/ qc_trim/ fastp.json -o multiqc_report/Open multiqc_report/multiqc_report.html. The General Statistics table at the top gives one row per file (reads, % duplicates, % GC, length). In the FastQC section, the trimmed files should show better per-base quality and almost no adapter content. Ready for alignment means: per base quality passes, adapter content is near zero, and you lost only a small share of reads.
Quick check: your raw reads show 18% adapter content at the last position. Roughly what share of your fragments are shorter than the read length?
About 18%. Adapter content at the final position equals the share of inserts shorter than the read length, because every such read runs into the adapter.
▸ Go deeper: getting more out of MultiQC
- MultiQC writes a
multiqc_data/folder next to the report.multiqc_general_stats.txtholds the General Statistics table as plain text, ready for R, Python or your supplementary tables. - Use
--filenameto name the report per project, and-fto overwrite an old one. - If sample names look messy (for example
sample_R1_fastqc), MultiQC’s sample-name cleaning and its config file let you rename samples or group read 1 and read 2. - Run MultiQC again after alignment and quantification: it also reads STAR, HISAT2, samtools, Picard, featureCounts and many more, giving one report for the whole project.
🔬 Researcher notes: set your acceptance criteria before you look
Decide what makes a sample pass before you see the results, then apply the rule to every sample. That avoids quietly keeping a sample you like and dropping one you don’t. A practical rule set for a sequencing batch:
- Read count within the range the facility promised (and similar across samples)
- Adapter content near 0% after trimming; less than about 5% of reads removed by trimming
- Q30 fraction typical of the instrument for every sample
- GC content within a few percent of the other samples of the same type
- Flag any sample that is an outlier on two or more metrics, and look at its full report before deciding
Step 6: Process many samples with one loop
Real projects have many samples. Put the raw files in raw/ named SAMPLE_R1.fastq.gz and SAMPLE_R2.fastq.gz, then run:
mkdir -p trimmed reports qc_raw qc_trim
fastqc raw/*.fastq.gz -o qc_raw/ -t 4
for r1 in raw/*_R1.fastq.gz; do
s=$(basename "$r1" _R1.fastq.gz)
fastp -i raw/${s}_R1.fastq.gz -I raw/${s}_R2.fastq.gz \
-o trimmed/${s}_R1.fastq.gz -O trimmed/${s}_R2.fastq.gz \
--detect_adapter_for_pe --cut_tail --cut_mean_quality 20 \
--length_required 36 -w 4 \
-h reports/${s}.fastp.html -j reports/${s}.fastp.json
done
fastqc trimmed/*.fastq.gz -o qc_trim/ -t 4
multiqc qc_raw/ qc_trim/ reports/ -o multiqc_report/▸ Go deeper: running samples in parallel
The loop handles one sample at a time. On a computer with many cores, run several samples at once with GNU parallel (install with conda install -c conda-forge parallel):
ls raw/*_R1.fastq.gz | sed 's#raw/##; s#_R1.fastq.gz##' | \
parallel -j 4 'fastp -i raw/{}_R1.fastq.gz -I raw/{}_R2.fastq.gz \
-o trimmed/{}_R1.fastq.gz -O trimmed/{}_R2.fastq.gz \
--detect_adapter_for_pe --cut_tail --cut_mean_quality 20 \
--length_required 36 -w 2 \
-h reports/{}.fastp.html -j reports/{}.fastp.json'Here 4 samples run together, each with 2 threads (8 cores in total). Match -j × -w to the cores you have.
🔬 Researcher notes: from scripts to workflows
For projects with dozens of samples, a workflow manager such as Snakemake or Nextflow re-runs only what changed, restarts after failures, runs on clusters and the cloud, and records exactly what was done. Established community pipelines (for example nf-core) already contain FastQC, trimming and MultiQC steps, so for standard data types it is often better to run a published pipeline and cite it than to maintain your own scripts. Learn to build them in our Reproducible Pipelines with Nextflow, Snakemake, Docker and HPC course.
Step 7 (optional): Check for contamination
FastQC tells you whether reads are well sequenced, not what organism they come from. Run a contamination check when GC content shows a second peak, when overrepresented sequences have no known source, or when the mapping rate later turns out lower than expected.
- FastQ Screen maps a sample of reads to several genomes (human, mouse, E. coli, PhiX, adapters, rRNA and so on) and reports the percentage hitting each one.
- Kraken2 classifies every read taxonomically using exact k-mer matches against a database, which is useful when you don’t know what the contaminant might be.
▸ Go deeper: running FastQ Screen and Kraken2
# FastQ Screen (download ready-made genome indexes once, about 10 GB)
conda install -c bioconda fastq-screen
fastq_screen --get_genomes
fastq_screen --conf FastQ_Screen_Genomes/fastq_screen.conf --outdir screen/ sample_R1.fastq.gz
# Kraken2 (needs a database, for example the standard or PlusPF database)
conda install -c bioconda kraken2
kraken2 --db k2_standard --threads 4 --paired --gzip-compressed \
--report sample.k2report --output /dev/null \
sample_R1.fastq.gz sample_R2.fastq.gzBoth tools produce reports that MultiQC can read, so contamination results appear next to your FastQC and fastp results.
🔬 Researcher notes: contamination, rRNA and sample identity
- FastQ Screen reports reads that map uniquely to one genome separately from those that map to several, which separates true contamination from conserved or low-complexity sequence (Wingett & Andrews 2018).
- Kraken2 uses minimizers and a compact hash table for speed and lower memory than the original Kraken (Wood et al. 2019). Classification is only as good as the database: reads from organisms that are missing from it stay unclassified or are assigned to relatives.
- For RNA-seq, check the rRNA fraction: a high value means poor rRNA depletion or poly-A selection, and trimming cannot fix it.
- Sample swaps are not visible in read QC. After alignment, check sex-chromosome coverage against the sample sheet and compare SNP genotypes between samples from the same individual.
Troubleshooting
| Problem | Likely cause | What to do |
|---|---|---|
fastqc: command not found | Environment not active | Run conda activate qc |
| FastQC runs out of memory | Very large files or many threads | Lower -t or add --memory 1024 (FastQC 0.12+) |
| Adapter content still high after fastp | Unusual adapter, single-end data | Give the adapter with --adapter_sequence (and --adapter_sequence_r2) |
| More than 20% of reads removed | Very poor run or settings too strict | Check the raw report; try --cut_mean_quality 15; ask the sequencing provider |
| Read counts in R1 and R2 differ after trimming | Files trimmed separately | Always trim both files in one fastp command |
| Two peaks in GC content | Contamination or adapter dimers | Check overrepresented sequences; screen with a contamination tool such as FastQ Screen or Kraken2 |
| MultiQC finds no reports | Wrong folder or old files | Point it to the folders with the _fastqc.zip and .json files |
| Long runs of G at read ends | Two-colour chemistry (NovaSeq/NextSeq): “no signal” is read as G | fastp trims poly-G automatically (force with -g); Cutadapt: --nextseq-trim=20 |
| Trimming removes almost everything, or nothing | Wrong quality encoding (Phred+64 data read as Phred+33, or the reverse) | Check Encoding in FastQC Basic Statistics; convert old data with seqtk seq -Q64 -V |
| Duplication FAIL on RNA-seq | Highly expressed genes; FastQC ignores the mate read | Expected. Assess duplication after alignment with paired-aware tools; do not deduplicate |