Quality Control Fundamentals

QC turns a set of files and tool reports into a documented decision about what can be analyzed, what needs investigation, and why. This lesson follows a paired-end short-read example, then contrasts it with single-cell checks. Thresholds must come from the assay and study context.

What you will produce

Work through the checks below and finish with a sample-level QC record: tool versions, selected metrics, a pass/review/exclude decision, its evidence, and a next action. The commands are templates, not a tested end-to-end pipeline.

• Define and interpret common quality-control (QC) metrics for sequencing experiments: per-base read quality, adapter content, mapping rate, read-pairing, duplication, coverage, contamination, and sample missingness.
• Differentiate technical failure from plausible biological signals and identify useful short-read, long-read, and single-cell QC views.
• Create a concise, reproducible QC record that lists tools and versions, key metrics, a pass/review decision, and the next action for each sample.
• Apply basic troubleshooting steps and document exclusions without deleting original raw data.

Prerequisites:

• Familiarity with command-line basics and file paths.
• Completed an introduction to biological data formats (FASTQ, BAM/CRAM, VCF, gene-count matrices).
• Thresholds depend on assay, organism, library protocol, instrument, and study design. This tutorial does not prescribe universal cutoffs.

Start with the question QC must answer

QC reports describe technical properties of the data; they do not establish whether a biological hypothesis is true. Poor base quality, contamination, low coverage, or duplication can affect downstream analyses. An unusual measurement can also reflect the sample or protocol. For example, a high mitochondrial RNA fraction needs interpretation alongside tissue, cell handling, and the distribution across cells and samples.

QC matters at multiple scales:

• Read-level: per-base quality, adapters, read length distribution.
• Alignment-level: mapping fraction, properly paired reads, insert-size distribution.
• Duplicate and coverage-level: PCR duplicates, uniformity of coverage for target capture.
• Sample-level: cross-sample contamination, identity swaps, missingness.
• Single-cell: cell complexity, mitochondrial fraction, doublets.

Paired-end short-read workflow

This example assumes two gzip-compressed FASTQ files for one sample, matching provider checksums, and (for alignment) a reference FASTA with compatible indexes. Create the top-level raw/, reference/, and logs/ inputs or directories as needed; the files shown are placeholders, not a supplied dataset. Run from the project directory after installing and testing the named tools. Replace sample names and paths. The snippets illustrate stages; they do not install tools, download references, or configure a compute environment.

Save each shell stage in a Bash script and run it with bash stage-name.sh from the project directory. Each script starts with set -euo pipefail, so a failed command or pipeline stops that stage. Do not continue to later stages after a failure, and do not treat an old output file as evidence that a new run succeeded.

Required inputs (example):

• raw/ sample_R1.fastq.gz, sample_R2.fastq.gz: primary raw FASTQ files.
• raw/checksums.txt: provider SHA-256 manifest, or raw/md5checksums.txt: provider MD5 manifest. Entries must name files relative to raw/, for example sample_R1.fastq.gz. Use the provider's actual names and algorithm; do not replace a provider manifest with locally computed hashes.
• reference/ genome.fa and its index files (if performing alignment).

1. Verify input files before analysis

What you do: confirm that each FASTQ matches the checksum provided by the data source. This step prevents downstream QC confusion caused by partial or corrupted downloads.

set -euo pipefail
# Input: provider manifest and FASTQ files beneath raw/
# Output: logs/checksum_validation.log
mkdir -p logs
if [[ -f raw/checksums.txt ]]; then
  (cd raw && sha256sum -c checksums.txt) 2>&1 | tee logs/checksum_validation.log
elif [[ -f raw/md5checksums.txt ]]; then
  (cd raw && md5sum -c md5checksums.txt) 2>&1 | tee logs/checksum_validation.log
else
  echo "Missing provider SHA-256 or MD5 manifest; stop and obtain it." >&2
  exit 1
fi

# Confirm each compressed FASTQ can be decompressed without error:
gzip -t raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz
echo "Provider checksums and gzip integrity checks passed."

Interpretation: all files should report OK and gzip -t should exit with code 0. If any check fails, re-download that specific file rather than editing the manifest. Save checksum_validation.log in logs/.

2. Compare read-level reports

What you do: run a per-file QC tool that reports per-base quality, adapter content, GC distribution, and overrepresented sequences. FastQC and fastp are common choices; run MultiQC to aggregate multiple FastQC outputs into a single report.

set -euo pipefail
# Input: raw/*.fastq.gz
# Outputs: qc/fastqc/ (HTML + raw), qc/multiqc_read_qc.html
mkdir -p qc/fastqc
fastqc -t 4 -o qc/fastqc raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz

# Aggregate into one report
cd qc
multiqc -n multiqc_read_qc.html -o . fastqc
cd ..

FastQC writes a report for each input file. MultiQC searches the supplied directory for recognized report files and compiles a combined report; review sample names and modules to ensure it found the intended inputs.

3. Trim only when a report supports it

Fastp can trim adapters and apply other filters, so this example is not a neutral “adapter-only” operation. Inspect the reports and defaults for your version, justify each change, and retain the original FASTQ. This command writes a new paired output and reports; it does not overwrite raw files.

set -euo pipefail
# Input: raw/*.fastq.gz
# Output: trimmed/sample_R1.trim.fastq.gz, trimmed/sample_R2.trim.fastq.gz
mkdir -p trimmed qc
fastp -i raw/sample_R1.fastq.gz -I raw/sample_R2.fastq.gz \
      -o trimmed/sample_R1.trim.fastq.gz -O trimmed/sample_R2.trim.fastq.gz \
      --detect_adapter_for_pe --trim_front1 0 --trim_tail1 0 \
      --thread 4 --html qc/fastp_report.html --json qc/fastp_report.json

Note: adapt parameters to your platform and chemistry. The HTML and JSON reports are written to qc/fastp_report.html and qc/fastp_report.json, respectively; keep these alongside the other QC reports.

4. Align, sort, and inspect the BAM

Align reads with your chosen aligner (bwa-mem2 for short-read DNA, STAR for RNA-seq, minimap2 for long reads). Here is a short-read DNA example using bwa-mem2 followed by samtools flagstat and Picard MarkDuplicates.

The indexing line writes companion files beside the FASTA and is needed only once for that exact reference and aligner version. Use a writable reference copy, or skip that line when matching bwa-mem2 indexes have already been prepared. A FASTA .fai alone is not a bwa-mem2 index.

set -euo pipefail
# Inputs: reference/genome.fa and trimmed/*.fastq.gz
# Outputs: indexed BAMs, flagstat/idxstats logs, duplicate metrics
mkdir -p aligned qc
bwa-mem2 index reference/genome.fa
bwa-mem2 mem -t 8 reference/genome.fa \
    trimmed/sample_R1.trim.fastq.gz trimmed/sample_R2.trim.fastq.gz \
    | samtools view -b -o aligned/sample.unsorted.bam -

samtools sort -@ 4 -o aligned/sample.bam aligned/sample.unsorted.bam
samtools index aligned/sample.bam

# Alignment QC
samtools flagstat aligned/sample.bam 2>&1 | tee qc/alignment_flagstat.txt
samtools idxstats aligned/sample.bam 2>&1 | tee qc/alignment_idxstats.txt

# Mark duplicates (Picard)
picard MarkDuplicates I=aligned/sample.bam O=aligned/sample.dedup.bam \
      M=qc/picard_markdup_metrics.txt CREATE_INDEX=true

flagstat summarizes alignment flag categories, including mapped reads, proper pairs, and records marked duplicate; interpret categories in light of the tool and input. idxstats reports counts by reference sequence after indexing. Picard MarkDuplicates writes metrics and marks duplicates in a new BAM. These outputs can flag questions for follow-up, but unusual contig counts alone do not prove contamination.

5. Check depth where the assay needs it

set -euo pipefail
# Input: aligned/sample.dedup.bam
# Output: qc/coverage_mean_depth.txt
mkdir -p qc
samtools depth -aa aligned/sample.dedup.bam | \
  awk '{sum+=$3; cnt++} END {if (cnt>0) print sum/cnt; else exit 1}' \
  > qc/coverage_mean_depth.txt

# Target breadth: fraction of each interval overlapped, not mean depth.
bedtools coverage -a targets/regions.bed -b aligned/sample.dedup.bam > qc/targets_coverage.tsv

The -aa denominator includes every reference position in the BAM header, including contigs with no mapped reads. With -a alone, wholly unused contigs can be omitted and the reported mean can be too high for a whole-reference claim. This whole-reference mean can be inappropriate for a targeted assay and expensive to calculate on a large genome. For a panel, define the intended target denominator and report per-target depth and breadth instead. The displayed bedtools coverage command reports overlap counts, covered bases, interval length, and covered fraction; it does not report mean depth. Record quality, duplicate, and overlapping-pair filters: samtools depth excludes marked duplicates by default but counts overlapping mates separately unless -s is requested.

6) Contamination and sample identity checks

Use lightweight taxonomic classifiers for microbial contamination (Kraken2) and fingerprinting tools (e.g., VerifyBamID2 or Peddy for human samples) for swaps/contamination. Keep results in logs/.

set -euo pipefail
# Kraken2 taxonomic profile
# Input: trimmed sample FASTQ, kraken2 database path
# Output: qc/kraken2_report.txt
mkdir -p qc
kraken2 --db /path/to/kraken2_db --threads 4 --report qc/kraken2_report.txt \
        --paired trimmed/sample_R1.trim.fastq.gz trimmed/sample_R2.trim.fastq.gz

# VerifyBamID2: human BAM, matching reference FASTA and SVD resource panel.
# Replace the absolute panel prefix below; its .UD, .mu and .bed files must exist.
# Bioconda executable name: verifybamid2; source builds may use VerifyBamID.
samtools faidx reference/genome.fa
mkdir -p qc/verifybamid2
(
  cd qc/verifybamid2
  verifybamid2 --BamFile ../../aligned/sample.dedup.bam \
    --Reference ../../reference/genome.fa \
    --SVDPrefix /absolute/path/to/build_matched_panel_prefix \
    > verifybamid2.log 2>&1
)
# Outputs: qc/verifybamid2/result.selfSM and result.Ancestry, plus the log.

VerifyBamID2 is not the legacy verifyBamID --bam --vcf interface. Obtain a compatible SVD resource panel and check assembly and contig naming against the BAM and FASTA before running it. Its default result files are written in the working directory, so use a separate output folder for every sample. Inspect FREEMIX in result.selfSM; a contamination estimate alone does not identify a swapped sample or establish sample identity. This is a documented command template, not an executed contamination analysis.

7) Single-cell specific QC

Start from a count matrix produced by your pipeline (for example, Cell Ranger or STARsolo). The R route below is the primary example; it reads a 10x Matrix Market directory, calculates per-cell features, counts and mitochondrial fraction, and draws distributions. The input is not included. Confirm feature naming and organism before using the mitochondrial pattern.

R workflow: Seurat QC metrics

# Input: a prepared 10x filtered_feature_bc_matrix/ directory
# Output: a Seurat object in memory and qc/seurat_qc_violin.png
library(Seurat)
library(ggplot2)
dir.create("qc", showWarnings = FALSE)
counts <- Read10X("data/filtered_feature_bc_matrix/")
seurat_obj <- CreateSeuratObject(counts = counts)
# Human gene-symbol convention only; adapt to species/annotation or omit.
mt_features <- grep("^MT-", rownames(seurat_obj), value = TRUE)
if (!length(mt_features)) stop("No features match ^MT-; adapt this pattern or omit the mitochondrial metric.")
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, features = mt_features)
qc_metrics <- seurat_obj[[]][, c("nFeature_RNA", "nCount_RNA", "percent.mt"), drop = FALSE]
print(summary(qc_metrics))
p <- VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
             ncol = 3, pt.size = 0.1)
ggsave("qc/seurat_qc_violin.png", plot = p, width = 9, height = 4)

Python alternative: Scanpy QC metrics

# Input: the same prepared 10x filtered_feature_bc_matrix/ directory
# Output: qc/scanpy_qc_violin.png and printed per-cell summaries
import os
import scanpy as sc
import matplotlib.pyplot as plt

os.makedirs("qc", exist_ok=True)
adata = sc.read_10x_mtx("data/filtered_feature_bc_matrix/")
adata.var_names_make_unique()
# Human gene-symbol convention only; adapt to species/annotation or skip.
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)
qc_columns = ["n_genes_by_counts", "total_counts", "pct_counts_mt"]
print(adata.obs[qc_columns].describe())
sc.pl.violin(adata, qc_columns, jitter=0.25, multi_panel=True, show=False)
plt.savefig("qc/scanpy_qc_violin.png", bbox_inches="tight")
plt.close()

Inspect distributions of genes and counts per cell, mitochondrial fraction, and relationships between metrics across samples. The mitochondrial pattern assumes human-style gene symbols; adapt it for the organism and annotation or omit it when inappropriate. Use plots to evaluate candidate filters; do not copy thresholds from another study without justification. These snippets calculate metrics and visualize them; they do not choose or apply a universal cell-filtering threshold.

Read the reports as comparisons, not verdicts

Below are common metrics, what they mean, and important caveats.

Per-base quality and adapter content

What to look for: consistently high Phred scores across reads; a low-quality tail at read ends suggests trimming. Adapter content above background indicates library preparation or sequencing artifacts.

Meaning and limits: high average quality reduces base-calling error, but does not guarantee correct alignment if the reference is distant. Adapter presence often indicates short inserts; confirm insert-size distribution before trimming aggressively.

For long reads, inspect read-length and per-read quality distributions (for example, with the platform's read-QC report) alongside mapped-read yield and coverage breadth/depth. A long-read length histogram can reveal short-fragment tails or a mixture of molecule sizes, but neither read length nor a single quality summary establishes assembly or variant-call accuracy. Compare samples prepared and sequenced with compatible chemistry, and consult platform-specific guidance before setting filters; this page does not provide a long-read processing recipe.

Mapping rate and properly paired reads

What to look for: mapping fraction and properly paired percentage in the context of the assay and comparable samples. A low mapping fraction or unexpected contig counts can prompt checks of reference build, library type, and possible contamination.

Meaning and limits: a low mapping rate can be caused by a wrong reference build, contamination, library type mismatch (e.g., RNA reads to DNA reference), or very divergent genomes. A high mapping rate alone does not guarantee correct variant calls; also check the mapping quality (MAPQ) distribution and variant-level metrics.

Duplicate fraction and library complexity

What to look for: high duplicate fraction indicates low complexity or over-amplification. Duplicate fractions depend on library preparation, sequencing depth, and assay; compare like samples and consider library complexity and UMI design.

Meaning and limits: Marked duplicates inflate apparent coverage. Removing duplicates is appropriate for many DNA analyses but may be inappropriate for some RNA-seq quantification or UMIs (unique molecular identifiers); treat UMI-aware processing differently.

Coverage uniformity and target dropouts

What to look for: per-target coverage reports that show consistent depth across the target set. Regions with zero or very low coverage are problematic and make variant calling unreliable at those loci.

Meaning and limits: average depth does not substitute for per-target depth. For clinical genes or actionable loci, document per-site coverage before making claims about negative findings.

Contamination and identity checks

What to look for: Kraken2 reporting a sizable fraction of non-target taxa; VerifyBamID reporting non-zero freemix suggestive of cross-sample contamination. Peddy and other tools can detect sample swaps using genotypes.

Meaning and limits: low-level microbial reads may reflect environmental DNA or index-hopping in multiplexed runs. Interpret contamination estimates against the method, controls, and downstream analysis. A single numeric estimate is not a universal acceptability threshold; document the evidence and its limits.

When a sample looks different

• Symptom: Unexpectedly low mapping rate
  Likely causes: Common causes include a wrong reference build, contamination, an unexpected sample type (e.g., a microbial sample aligned to a human reference), a degraded or low-complexity library, adapter contamination, or very short reads.
  Fixes: confirm reference fasta and its index; run Kraken2 to detect major taxa; inspect FastQC for adapters/short inserts; re-run alignment with relaxed settings if reads are from a divergent strain; consider re-sequencing.
• Symptom: Unexpectedly high duplicate fraction
  Likely causes: low input DNA, over-amplification, targeted enrichment, or issues with UMI handling.
  Fixes: check library prep notes, use UMI-aware duplicate handling if UMIs were used, evaluate whether duplicates should be marked or not for your intended analysis; for low-input samples consider re-preparing with higher input if possible.
• Symptom: High mitochondrial RNA fraction in scRNA-seq cells
  Likely causes: abundant dead/dying cells, dissociation-induced stress, or tissue-specific biology (e.g., muscle).
  Fixes: inspect distribution across samples/batches; if confined to a subset consider excluding those cells but document criteria and numbers excluded; consider alternative dissociation protocols in future experiments.
• Symptom: Unexpected taxa dominate Kraken2 report
  Likely causes: environmental contamination, index bleed, or mislabeled sample.
  Fixes: check lab records and sample sheet for potential swaps; examine negative controls; if the contamination compromises analysis, remove contaminated libraries or re-sequence.
• Symptom: Large differences in QC metrics between batches
  Likely causes: library preparation batch effects, different sequencer models/chemistries, operator differences.
  Fixes: include batch covariates in downstream models, consider re-processing to harmonize protocols, or re-sequence critical samples if possible.

Record decisions and protect source data

Tools and databases change. Always record tool names and exact versions (e.g., fastp v0.23.2, bwa-mem2 v2.2.1, samtools v1.17). Reference genomes and annotation versions matter: chromosome names (chr1 vs 1), contig order, and the presence of decoy or alternate contigs affect mapping statistics. When using taxonomic classifiers (Kraken2), record the database build date and the specific database used. Contamination detection is only as good as the database.

Responsible-use notes:

• Do not delete raw files. Keep originals in raw/ and store logs alongside them.
• Document all QC decisions: who made them, when, and why. If samples are excluded, record the exact rationale and provide the outputs that led to exclusion.
• Do not apply rigid thresholds blindly. Use data visualization to understand distributions and make context-aware decisions.
• Privacy: human sequencing data are sensitive. Ensure appropriate access controls when storing and sharing QC outputs that may contain sample identifiers.

Make the QC record actionable

Before interpreting a flag: compare samples processed with the same assay and reference, check the report’s sample identity, and review controls and batch structure. A flagged metric is a prompt to investigate, not an automatic exclusion rule.

• Store provider checksum manifest and logs in logs/ and confirm all files are OK.
• Run FastQC (or fastp) on raw FASTQ and aggregate with MultiQC; save the HTML report.
• Decide whether trimming is justified; if so, produce trimmed FASTQ in a separate directory and save trimming reports.
• Align reads, index BAMs, run flagstat and Picard MarkDuplicates; store outputs under qc/.
• Compute coverage metrics and per-target coverage if applicable.
• Run contamination/identity checks and save reports.
• For single-cell data, compute per-cell QC metrics and visualize distributions before filtering.
• Create a QC summary file (e.g., qc/qc_summary.tsv) listing sample, tool versions, key metrics, decision (pass/review/exclude), and next action.

Use one row per sample and write NA when a metric was not measured; do not treat missing evidence as a pass. The following row is illustrative, not a universal cutoff or a biological benchmark. Its review decision follows from the stated study-specific comparison: the sample's mapped fraction is below its matched batch range, so reference compatibility and read-level reports should be checked before downstream interpretation.

sample	tool_versions	metrics	decision	evidence	next_action
SYNTH_A	fastp 0.23.2; samtools 1.17	reads=1,000,000; mapped=62%; proper_pairs=58%; mean_depth=18x	review	mapped fraction below matched-batch range (illustrative); controls pending	verify reference build and inspect FastQC/MultiQC reports; reassess with controls

To make the record reproducible, place it at qc/qc_summary.tsv, retain the corresponding reports in qc/, and record exact installed versions rather than copying the illustrative version labels. A completed row documents a decision and its evidence; it does not by itself validate a sample or establish a universal pass threshold.

What this run produced

This is the actual Seurat QC plot saved by the R example, using the official 10x Genomics PBMC3k filtered matrix. It is a public single-cell example, not the result of the FASTQ, alignment or bulk-QC templates on this page.

Run details and limits

The exact R block ran with R 4.3.3, Seurat 5.0.1, SeuratObject 5.0.1, Matrix 1.6.5, and ggplot2 3.4.4 on Linux. It read 32,738 features and all 2,700 supplied cells without additional sampling or filtering.

The documented human ^MT- feature pattern was found, the three metric summaries printed, and the PNG was saved. This checks only the prepared-input R component. Raw-read QC, bulk assays, other species, and the Python alternative were not executed by this run. Inspect the execution evidence.

Official tool documentation

• FastQC: Babraham Bioinformatics
• MultiQC documentation
• fastp documentation
• SAMtools documentation
• SAMtools depth: denominator and filtering options
• VerifyBamID2: reference-panel prerequisites and command interface
• VerifyBamID2 wrapper: executable name and default result files
• GNU SHA checksum documentation
• Scanpy QC metrics API
• Seurat PercentageFeatureSet reference
• Picard tools documentation
• Kraken2: ultrafast metagenomic sequence classification
• Scanpy documentation (single-cell analysis in Python)