FASTQ QC Pipeline

Build a complete quality-control pipeline for FASTQ files. You will learn how to read sequencing data, compute quality statistics, filter reads, count k-mers, and generate a summary report.

Prerequisites

Complete the Hello Genomics tutorial first. You will also need a sample FASTQ file. BioLang ships with sample data:

# Sample data is included in the repository
ls examples/sample-data/reads.fq

# Or run the quickstart to verify everything works
bl run examples/quickstart.bl

The sample FASTQ file at examples/sample-data/reads.fq contains 8 reads with varying quality scores. For this tutorial, we will use it directly or you can substitute your own FASTQ file.

Step 1 — Reading a FASTQ File

Use read_fastq() to load a FASTQ file into a reusable table, or fastq() for a lazy stream (better for large files). Each record contains an id, sequence, and quality scores.

# requires: examples/sample-data/reads.fq in working directory
# qc.bl — FASTQ quality control pipeline

# read_fastq returns a Table (reusable)
let reads = read_fastq("examples/sample-data/reads.fq")

# Check how many reads we have
println("Total reads:", len(reads))

# Look at the first read
let first = reads[0]
println(f"ID:       {first.id}")
println(f"Sequence: {first.seq}")
println(f"Quality:  {first.quality}")
println(f"Length:   {first.length}")

bl run qc.bl
# Total reads: 100
# Header:   @SRR001/1
# Sequence: ATCGATCGATCG...
# Quality:  IIIIIIHHGGFF...
# Length:   150

Step 2 — Computing Read Statistics

Use read_stats() to compute summary statistics across all reads, including quality scores, GC content, and length distribution.

# Compute read-level statistics using read_stats()
# read_stats returns a record with count, mean_length,
# mean_quality, gc_content, q20_pct, q30_pct, etc.
let stats = read_stats("examples/sample-data/reads.fq")

println(f"Total reads:  {stats.total_reads}")
println(f"Mean quality: {round(stats.quality.mean, 2)}")
println(f"Mean length:  {round(stats.length.mean, 1)}")
println(f"Mean GC:      {round(stats.gc_content, 4)}")

# You can also look at individual reads
let first = reads[0]
println(f"Read 1 GC: {gc_content(first.seq)}")
println(f"Read 1 length: {first.length}")

Step 3 — Exploring Per-Read Metrics

Beyond aggregate statistics, you can compute per-read metrics like GC content using map() and standard stats functions.

# read_stats() computes comprehensive statistics in one call
let stats = read_stats("examples/sample-data/reads.fq")

println("=== Read Statistics ===")
println(f"Total reads:  {stats.total_reads}")
println(f"Mean length:  {round(stats.length.mean, 1)}")
println(f"Mean quality: {round(stats.quality.mean, 2)}")
println(f"Mean GC:      {round(stats.gc_content, 4)}")
println(f"Q20 %%:        {round(stats.quality.q20_pct, 1)}")
println(f"Q30 %%:        {round(stats.quality.q30_pct, 1)}")

# You can also compute GC for individual sequences
let gc_values = map(reads, |r| gc_content(r.seq))
println(f"GC range: {min(gc_values)} - {max(gc_values)}")

Step 4 — Filtering Reads by Quality

Low-quality reads introduce errors in downstream analysis. Let us filter them out.

# filter_reads(input_path, output_path, options)
let result = filter_reads(
  "data/reads.fastq",
  "data/reads.filtered.fastq",
  {
    min_quality: 20,
    min_len: 50,
    max_n: 5,
    complexity: 0.3
  }
)

println(f"Passed filters: {result.passed} of {result.total} reads")
let clean = read_fastq("data/reads.filtered.fastq")

Step 5 — Trimming Low-Quality Ends

Instead of discarding entire reads, we can trim low-quality bases from the 3' end.

# trim_quality(input_path, output_path, threshold) trims the 3' end
let trim_result = trim_quality(
  "data/reads.fastq",
  "data/reads.trimmed.fastq",
  20
)
let trimmed_list = read_fastq(trim_result.output)
println(f"After trimming: {trim_result.passed_reads} reads")

# Compare stats before and after
let raw_stats = read_stats("data/reads.fastq")
let trim_stats = read_stats(trim_result.output)
println(f"Mean length before: {round(raw_stats.length.mean, 1)}")
println(f"Mean length after:  {round(trim_stats.length.mean, 1)}")

# Sliding-window trimming can also enforce a minimum length
let window_result = trim_reads(
  "data/reads.fastq",
  "data/reads.window-trimmed.fastq",
  {quality: 20, min_len: 30, window: 4, cut_front: true, cut_tail: true}
)
println(f"After window trimming: {window_result.passed_reads} reads")

Step 6 — K-mer Counting

K-mer profiles help detect contamination, adapter sequences, and biased composition. kmer_count accepts lists, tables, and streams directly — no need to extract sequences manually.

# Count k-mers across all trimmed reads
# kmer_count returns a Table sorted by count (descending)
let k = 5
let kmer_counts = kmer_count(trimmed_list, k)
println(f"\n=== Top {k}-mers ===")
kmer_counts |> head(10) |> println()

# You can also get k-mers for a single sequence
let first_kmers = kmers(trimmed_list[0].seq, k)
println(f"K-mers in first read: {len(first_kmers)}")

# Check for adapter contamination
let detected = detect_adapters("data/reads.fastq")
println("Detected adapters:")
println(detected)

Step 7 — Building the Complete Pipeline

Now let us combine everything into a single pipeline using the pipe operator.

# requires: examples/sample.fastq in working directory
# fastq_qc.bl — complete FASTQ QC pipeline

fn run_qc(input_path, output_path) {
  println(f"Reading {input_path}...")
  # Trim and filter in one file-to-file pass
  let trimming = trim_reads(input_path, output_path, {
    quality: 25,
    min_len: 50,
    window: 4,
    cut_front: true,
    cut_tail: true
  })
  let clean = read_fastq(output_path)

  # Compute stats on both raw and clean
  let raw_stats = read_stats(input_path)
  let clean_stats = read_stats(output_path)

  # Print comparison
  println("\n=== QC Report ===")
  println("                  Raw        Clean")
  println(f"Reads:        {raw_stats.total_reads}     {clean_stats.total_reads}")
  println(f"Mean length:  {round(raw_stats.length.mean, 1)}  {round(clean_stats.length.mean, 1)}")
  println(f"Mean quality: {round(raw_stats.quality.mean, 2)}     {round(clean_stats.quality.mean, 2)}")
  println(f"Mean GC:      {round(raw_stats.gc_content, 4)}   {round(clean_stats.gc_content, 4)}")
  println(f"Q30 %%:        {round(raw_stats.quality.q30_pct, 1)}      {round(clean_stats.quality.q30_pct, 1)}")

  println(f"\nWrote {trimming.passed_reads} clean reads to {output_path}")
}

# Run the pipeline
run_qc("examples/sample.fastq", "examples/sample_clean.fastq")

Step 8 — Generating a Summary Report

Build a comparison table of raw vs clean statistics and export it as CSV.

# requires: examples/sample.fastq and examples/sample_clean.fastq in working directory
# Generate a summary table comparing raw vs clean reads
let raw = read_fastq("examples/sample.fastq")
let clean = read_fastq("examples/sample_clean.fastq")

let raw_stats = read_stats("examples/sample.fastq")
let clean_stats = read_stats("examples/sample_clean.fastq")

# Build a comparison table
let labels = ["count", "mean_length", "mean_quality", "gc_content", "q20_pct", "q30_pct"]
let raw_vals = [raw_stats.total_reads, round(raw_stats.length.mean, 1),
                round(raw_stats.quality.mean, 2), round(raw_stats.gc_content, 4),
                round(raw_stats.quality.q20_pct, 1), round(raw_stats.quality.q30_pct, 1)]
let clean_vals = [clean_stats.total_reads, round(clean_stats.length.mean, 1),
                  round(clean_stats.quality.mean, 2), round(clean_stats.gc_content, 4),
                  round(clean_stats.quality.q20_pct, 1), round(clean_stats.quality.q30_pct, 1)]

let summary = to_table(map(
  [0, 1, 2, 3, 4, 5],
  |i| {metric: labels[i], raw: raw_vals[i], clean: clean_vals[i]}
))

println(summary)
write_csv(summary, "examples/qc_summary.csv")
println("Summary saved to examples/qc_summary.csv")

bl run fastq_qc.bl
# Reading examples/sample.fastq...
# === QC Report ===
# ...
# Summary saved to examples/qc_summary.csv