genome-pipeline

How it is built

Architecture overview

A genome pipeline turns the raw output of a sequencer, hundreds of millions of short text reads of a person's DNA, into a list of the places where that person differs from a reference genome. It does this in four stages: align each read to the reference, organise the alignments (sort, index, mark duplicates), call variants by comparing the reads piled up at each locus, and write the result as a compact, indexed file format. The idea that makes it tractable is seed and extend: never compare a read to the whole genome, find a few exact short matches through an index, and run an expensive alignment only around them. The same separation applies downstream, where the expensive probabilistic model runs only on regions where the reads disagree with the reference. BWA and minimap2 are the aligners, samtools and htslib the file layer, and GATK the best-known caller.

Parts:

• Reference index — find where a short string occurs in a multi-gigabase reference — fm-index (no page yet), built on Burrows-Wheeler transform, Suffix array, Rank and select bitvectors
• Minimizer index — the alternative seeding index for long or noisy reads — minimizer-kmer-index (no page yet), built on rolling-hash (no page yet), Swiss table, Fast hash functions
• Aligner — extends seeds into scored base-level alignments — Smith-Waterman alignment, SIMD programming, edit-distance (no page yet)
• Alignment files — SAM/BAM/CRAM and VCF records on disk, with random access — seekable-compressed-stream (no page yet), Huffman coding, run-length-encoding (no page yet), Interval tree
• Sort and merge — put hundreds of GB of alignments in coordinate order — External merge sort, k-way-merge (no page yet), pdqsort, Radix sort, Buffer pools (zero-alloc I/O), fsync and direct I/O
• Local assembly — rebuild candidate haplotypes where reads disagree with the reference — de-bruijn-graph-assembly (no page yet), built on Bloom filter, BFS and DFS
• Variant likelihoods — score each haplotype against each read, then each genotype per site — pair-hmm (no page yet), bayesian-genotype-likelihoods (no page yet), Floating-point error
• Variant files — normalised VCF records, joint across samples — vcf-variant-normalization (no page yet), Sparse matrices (CSR)
• Parallel runtime — threads, scatter-gather over genome intervals — Memory ordering

 FASTQ reads ─► seed: index lookup ─► chain / extend ─► SAM records (one per read)
                (FM-index | minimizers)  (Smith-Waterman, SIMD)       │
                                                                       ▼
        sort by coordinate ◄── BGZF-compressed BAM chunks ◄── write (threads, buffers)
        (external merge sort)                  │
              │                                ▼
              ├─► mark duplicates ─► analysis-ready BAM + index (.bai/.csi)
              ▼
   per interval (scatter): pileup ─► active regions ─► local assembly ─► pair-HMM ─► genotype likelihoods
              │
              ▼
   GVCF per sample ─► joint genotyping across samples ─► VCF (normalised, indexed) ─► filter / annotate

A whole genome is processed by scatter: the reference is cut into intervals, each interval is a job, and the results are gathered. The design has two throughput regimes: the aligner is CPU and memory-latency bound on one machine, the rest is I/O and sorting bound across files.

How the parts fit together

Reference index: two designs

The aligner needs an index of the reference that answers "where does this 20-mer occur?". Two designs compete, and they are alternatives:

• FM-index (fm-index, no page yet): a Burrows-Wheeler transform of the reference plus rank structures supports backward search, counting and locating a substring in compressed space, and the Suffix array (sampled) turns a match range into positions. BWA builds the FM-index of the reference with its index command (source: BWA README). It finds exact matches of any length, which suits short, accurate reads and keeps the index small for a human-sized genome.
• Minimizer index (minimizer-kmer-index, no page yet): from each window of w consecutive k-mers keep only the smallest by hash, store the kept k-mers in a hash table (Swiss table with a fast hash; rolling-hash, no page yet, is the cheap way to hash overlapping k-mers). Minimap2 uses this: it extracts (k,w)-minimizers of the reference and indexes them in a hash table, ignores the most frequent minimizers (the -f option, 2e-4 by default), and reads reference bases in batches of 4G by default (source: minimap2 README, algorithm overview). A table lookup per minimizer is faster than a backward search, and sampling makes seeds fewer, so it scales to long, noisy reads (PacBio, Nanopore) where exact long matches are rare. The BWA README states that BWA-MEM has been replaced by minimap2 for long reads. The contract with the aligner is the same in both cases: a query k-mer in, a list of reference positions out, with repeat-heavy seeds suppressed.

Aligner

Seeds are not alignments. Both aligners collect seeds, then chain them: minimap2 sorts seeds by reference position and chains them with dynamic programming, and each chain is a candidate mapping (source: minimap2 README). Only then does the expensive step run, base-level alignment of the gaps between seeds with Smith-Waterman alignment (affine gaps, O(nm)), usually banded and vectorised with SIMD programming: many cells or many reads per instruction. Bit-parallel edit distance (edit-distance, no page yet) serves where only the number of differences matters. Minimap2 splits chains where the accumulated score drops by -z (default 400) and filters seeds likely to make both a long insertion and a long deletion. The output is one SAM record per alignment (source: SAM spec) with the position, a mapping quality, a CIGAR string for the base-level path and flags for pairing, strand and primary or secondary status.

Alignment files

A SAM line has eleven mandatory fields (name, flag, reference, position, mapping quality, CIGAR, mate fields, sequence, quality). The binary form BAM is compressed in BGZF, a gzip variant made of independent blocks, each at most 64 KiB before and after compression (source: SAM spec, section 4.1). Independent blocks give random access: an index stores virtual offsets, an unsigned 64-bit value of block start shifted left 16 bits plus the offset inside the block (coffset<<16 | uoffset; source: SAM spec). This is the seekable-compressed-stream block (no page yet); Huffman coding is the entropy stage of the DEFLATE inside each block. CRAM is a further format for the same data. The htslib README describes the C library that reads SAM, CRAM and VCF and is the core of samtools and bcftools; it implements a generalised index with extension .csi and ships the tabix indexer and the bgzip utility (source: htslib README). Queries "all reads overlapping chr7:1,000,000-1,001,000" go through a binning or tree index; Interval tree is the in-memory version of the same idea, used for annotation and interval overlap. The RLE and bit-packing blocks fit qualities, which are long runs of similar small integers (run-length-encoding, no page yet).

Sort and merge

Callers need reads in coordinate order so that all reads over a locus are adjacent; the SAM header records this as SO:coordinate (source: SAM spec). Sorting a BAM bigger than memory is External merge sort: sort memory-sized chunks (pdqsort, or Radix sort on the integer coordinate key), spill each as a compressed temporary file, then k-way merge (k-way-merge, no page yet) them. Buffer pools keep the hot path from allocating, and the I/O rules in fsync and direct I/O decide how temporary files are written. Radix beats comparison sort when the key is a small integer; pdqsort is the safe default for composite keys such as (reference, position, strand). Duplicate marking works on the sorted file because duplicate fragments from the same molecule land at the same coordinates (flag 0x400 marks them; source: SAM spec).

Local assembly and likelihoods

Single-nucleotide differences are easy to find by counting. Insertions, deletions and clustered changes are not, because aligners place them inconsistently. A caller like GATK's HaplotypeCaller therefore finds active regions where reads disagree with the reference, reassembles the reads of that region on a graph of k-mers (de-bruijn-graph-assembly, no page yet: a Bloom filter or hash table counts k-mers, BFS and DFS walk paths; tips and bubbles are pruned), and extracts candidate haplotypes. Each read is scored against each haplotype with a pair hidden Markov model (pair-hmm, no page yet), a forward algorithm over match, insertion and deletion states vectorised with SIMD. Those read-haplotype likelihoods feed bayesian-genotype-likelihoods (no page yet): combined with a prior they give a likelihood for every genotype at a site. Products of many small probabilities underflow, so the arithmetic runs in log space or with scaling, which is the point of Floating-point error. The contract between these stages is a list of candidate alleles per site and a table of log-likelihoods (read by haplotype).

Variant files and joint calling

Output is VCF. Because the same indel can be written in several equivalent ways, vcf-variant-normalization (no page yet) left-aligns and trims indels so equal variants have one form, which is required before comparing or merging callsets. GATK's germline workflow runs HaplotypeCaller in GVCF mode on a single sample, scattered over an interval list, and produces one GVCF per sample; several GVCFs go to GenomicsDBImport and then to joint genotyping, producing a multi-sample VCF (source: gatk-workflows gatk4-germline-snps-indels README). The benefit of GVCFs is that adding a sample does not require re-reading every other sample's reads. A site-by-sample table is mostly "no variant", so Sparse matrices (CSR) are the natural layout (GenomicsDB's internal format may differ).

Parallel runtime

Threads in an aligner share the index read-only, so the cost is memory latency, not locking; the places that need care are the work queue and the ordered output writer, where Memory ordering decides how a producer hands a batch to a consumer. Across machines the unit is the genome interval and the sample, as in the scatter above.

Data and control flow

One sample, reads to variants:

1. FASTQ reads in; an index is already built once for the reference (bwa index, or minimap2 on the fly or a saved .mmi).
2. Each worker takes a batch of reads, looks up seeds, chains and extends; minimap2 reads query bases in batches of 200M by default (the -K option; source: minimap2 README).
3. Alignments are written as SAM records, compressed to BGZF blocks, and spilled to disk in sorted chunks.
4. Chunks are merged into a coordinate-sorted BAM; an index is written; duplicates are marked.
5. The genome is cut into intervals; per interval, the caller scans for active regions, assembles, runs the pair-HMM and writes GVCF records.
6. GVCFs are merged per interval across samples and jointly genotyped; the VCF is normalised and indexed. Time goes to the aligner's extension step and to the pair-HMM (both vectorised dynamic programming), and to I/O during sorting.

One lookup in the index (seed step):

1. Take a k-mer or minimizer from the read.
2. FM-index: backward search narrows a suffix-array interval character by character; minimizers: one hash lookup.
3. Convert the interval or bucket to reference positions; drop seeds from over-frequent k-mers.
4. Hand the positions to chaining.

Scaling

1. Alignment throughput is bound by memory latency of index lookups and the extension DP. Levers: more threads on one node, SIMD extension, cache-friendlier index layouts (BWA-MEM2 is a reimplementation whose output is meant to match BWA-MEM 0.7.17; source: its README), then more nodes with sample-level parallelism.
2. Sorting and file I/O are next: big temporary files, many small writes. Levers: larger sort buffers, more compression threads, fast local scratch, buffer pools, fewer intermediate files (stream aligner output into the sorter).
3. Variant calling parallelises across intervals nearly perfectly; the pair-HMM is the hot loop. Levers: SIMD or GPU/FPGA kernels for the HMM, tighter active regions, scatter into more intervals.
4. Cohort size: joint calling touches every sample at every site. The lever is the GVCF design plus a database import step, then splitting by interval. Order of: a human genome is about 3 billion bases, and 30x coverage means on the order of 100 billion sequenced bases, so every constant factor matters.

Failure modes

Real implementations

The alternative: a reference-free route assembles reads into contigs (de Bruijn graph or overlap graph) and compares assemblies, which finds large structural differences and novel sequence but costs more compute and memory than mapping. Pick mapping-plus-calling for resequencing one individual against a good reference; assemble when no reference exists.

Sources

• H. Li, Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM, arXiv 2013. paper
• H. Li, Minimap2: pairwise alignment for nucleotide sequences, Bioinformatics 2018. doi; algorithm overview in the minimap2 README.
• BWA README: github.com/lh3/bwa.
• SAM/BAM specification (BGZF, virtual offsets, sort order, flags): SAMv1.pdf; VCF: VCFv4.3.pdf.
• htslib and samtools: htslib README, htslib.org.
• GATK: repository and the germline workflow gatk4-germline-snps-indels.